SlffiffiB HHBHHHi |H|l|BHHBP|BRnOBBH|D]|ili^U awli BTnwnWBffltiTTffflawr |RHRHK8SalBlilH llHi Hfll ISttB&lNulnX IrnrtlTBmlHmflBlltiTOfPifflTnrTTi FHnlrTBTiii f' n * l *f LIuMuiJiill! »■» mmi wnmrw itHli BKffiKSs VHHBIBnn) ■ WHXBBBBH&m HHBH SpHi liiill phB HI HnffiBi IBBnfflffl MflHBHfiflHHnB9 M I ' IflUlffiB&YMltz ■BMHa Bfl Sfinn HhBE .«fe^>* K119IBIHI SHgsgl Hi H 111 ISIS LIBRARY OF THE UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN cop. a Digitized by the Internet Archive in 2013 http://archive.org/details/montecarlosimula519wang /^JIl IUCDCS-R-72-519 y?i^C4i MONTE CARLO SIMULATION OF NONLINEAR RADIATION INDUCED PLASMAS by Benjamin Shaw-hu Wang May, 1972 IH6: imm* 1) at the expense of more processing time. The reason for this is based on the fact that the system being simulated is time independent. We have introduced the time variable artificially so that ergodicity holds, namely ensemble averaging is the same as time averaging. 20 2. Tracing Procedure The simulation is initialized by introducing N particles with energy E. [sampled from an initial guess for the number density distribution vs. energy], with spatial position z. (uniformally distributed throughout the spatial region), with direction u. (u=cosfi; 9: angle of flight with z axis) which is uniformly distributed between -1 and +1 and with weight wt . . Volume particles are emitted uniformly throughout the medium and isotropically in direction; namely, at t=t : o E° 1 = E. in Z° 1 = r • D i 1 = 2r.-l l / (4a) (4b) (4c) where r.: random number, uniformly distributed between 0, 1, or x e (0, 1) d: slab width or cylinder diameter wt.,E. : weight and particle energy determined from initial distribution. Once the particle data are determined, the geometry routine is entered. This includes the determination of the distance to boundary d„ , the distance to collision d _ , and J B' col' the distance to census time d (clock distance equivalent cen ^ to At). These distances are defined by: 21 d B = [(j+l)AZ - Z]/u if u > (5) = (Z - jAZ)/ |u| if u < = DLB (a large number DLB » AZ) if u^r d CQl = A(E) |lnr| = |lnr| / J t (E) (6) cen D • DI = V • At* DI (7) avg avg where r 3 It D avg random number uniformly distributed zone number spatial interval size the total cross-section evaluated at the particles' energy the average value of D ^ of N particles of the v col * previous cycle. It is approximately one mean free path corresponding to E^ tt/ ^: the average energy of the cycle. avg DI : input parameter for controlling the cycle time interval. The next step depends on which of these distances is the smallest. a . The System Without Electric Field (E =0 ) i) If the smallest distance is d : the particle cen ^ is advanced in position and time to census time. z 1 = z + u.d u' = u E« = E cen (8a) (8b) (8c) The particle data is stored, the energy is tallied, and 22 the density counter is incremented. Then we go on to process the next particle. ii) If the smallest distance is d the particle is advanced to the boundary such that: z' = z + (d + 0.001«Az) . u j = INT (z»/Az) + 1 d T d' , = d , col col B d' = d - d n cen cen B (9a) (9b) (9c) (9d) where j is the zone number. If j>JZ or jD r z'^0 then the particle is considered to have escaped from the system. If DCEN > DCOL, the particle is advanced to collision point. At' = At - DCOL-At/DCEN DCEN = DCEN - DCOL Z' = z + Az-DCOL/DCEN E« = E + Az.E -DCOL/DCEN (22) then the specific collision type must be determined. 27 c. Types of Collisions and Cross Sections For the elastic collisions, the particle scatters off with angle of deflection 6 and energy loss &E, the follow- ing formula is used (36). AE = 4r < E ~ E J (1-cosG ) (23) M n s where 9 : scattering angle (angle of deflection) E: particle energy E : target gas atom energy m: electron mass M: gas atom mass E is selected according to a Maxwellian distribution, and n * ' the AE so formed may have either positive or negative values, indicating the loss or gain of energy respectively For calculational simplicity, instead of selecting the deflection from the elaborate distributions, we use 7T/4 or tt/2 deflection approximations, namely 9 = X/4 or V2, then cosB = or 0.70711. Another reasonable alternative s is to assume that the direction after collision is equally probably in all directions (isotropic scattering). We select this direction at random with uniform distribution, and based on Eq. 23 calculate cos9 and AE, respectively, where cos9 = 2r.-l (24) l 6=9 --6 . s 1 28 where 9 . : incident angle 1 ' 9 : outgoing angle after collision For the three kinds of inelastic collision cross sections we use approximate formula based on experimental data as correlated by Itoh and Musha (2). The formulas are expressed in the form: B(E-E. ) £( E ) = 1 j + C (25) A+(E-E ± ) Z where the constants A, B, C, and E. depend on the kind of collision and scattering angle. For the elastic cross section in helium, we use a single approximate formula for whole effective energy range based on Heylen and Lewis's paper (10), namely: T el (E) = 26exp(-0.04E) (cm -1 ) , (26) For the recombination cross section for helium (35), we assume the simple form: 1rcb (e) =^ 1, the system has more particles introduced 32 into than eliminated, but if BF < 1, the opposite is true. The over-correction factor is proportional to (BF) t i.e., to the kth power of the balance factor. The value of k. initially chosen was k=2, but an optimum value can be obtained by trial and error. 33 CHAPTER III VARIANCE REDUCTION TECHNIQUES A. General Background The Monte Carlo method (or method of statistical trials) consists of solving various problems of computational mathematics by the construction of some random process for each such problem, with the parameters of the process set equal to the required quantities of the problem. These quantities are then approximated by observations of the random process and the computation of its statistical characteristics, which are approximately equal to required parameters [from Schreider (15)]. Every Monte Carlo computation that leads to quantitative results may be regarded as estimating the value of a multiple integral. Suppose we have M random numbers in the computation system. The results will be a vector valued function R(^ 1 ,£?, ... £ M ) involving the sequence of random numbers £.. , £~, ... c M . This is an unbiased estimator of ... R(x, ,x~. ... x„) Jo Jo 1 2 M dx. ... dx . For the sake of simplicity, we take the one- dimensional integral as a standard example. Let r Jo f(x)dx (34) 2 f x 2 where f e L (0,1) or [f(x)] dx exists. / Define the relative efficiency of two Monte Carlo methods as follows. Let the methods call for n, and n ? units of computing time, respectively, and let the resulting 34 2 2 estimates of 9 have variances aC and 2 (35) Equation 35 is the product of two terms, the variance ratio n. CT^ 2 / 6^2 and the labor ratio "l/n 2 . If % , %, ?N are independent random numbers with uniform distribution between 0, 1 or *1 . e(0,l); then the quantities f.=f(^.) are the independent random variates with expectation 8 . Therefore N f4l f- N i=l X is an unbiased estimator of 6 , and its variance is (36) i r 1 (f(x)- J B) 2 dx= £ (37) Accordingly the quantity f~, which has been determined by observation of the random process, is approximately equal to the required quantity 8 (or in other words, it is an unbiased estimator of 8 ) with .a probability which can be made as close as required to unity if a large number of trials is practical. We refer to the estimator 7 as the crude Monte Carlo estimator of 8 (commonly referred to as Monte Carlo estimation) . Let us introduce another even less efficient method, namely hit or miss Monte Carlo. Suppose that < f (x) < 1 35 when < x < 1. We may draw the curve y = f(x) in the unit square < x, y < 1, and 9 is the proportion of the area of the square beneath the curve. Stating this formally, we write f(x) = g(x,y)dy (38) Jo where g(x,y) = if f(x) < y = 1 if f (x) > y '■:: . Then fl Jo JO The estimator n g(x,y) dxdy (39) 5-5? w'^n, <2i ) -5 1 = 1 (40) n* is the number of occasions where f( £ ~ . .. ) > £_. Y 2i-l — 7 2i In other words, we take n trials at random in the unit square, and count the proportion of them which lie below the curve y = f(x). Historically, a hit or miss method was first propounded in the explanation of Monte Carlo techniques. It is the easiest to understand but one of the least efficient techniques. The rate of convergence and the variance (statistical error) is in general proportional to /JN where N is the number of samples. It is obviously imprac- tical to gain a significant improvement of accuracy by merely increasing N. Therefore more efficient techniques for reducing the variance and/or increasing the efficiency 36 are very important. Monte Carlo experimentalists need wide experience and background in both the mathematics and the physics of the problem, and they have to exercise considerable ingenuity in distorting and modifying problems in the pursuit of variance reduction techniques. Statistical and inferential procedures are also important in order to extract the most reliable conclusions from the observational data. The basic requirement of the Monte Carlo calculations is to establish how random numbers may be used to sample a function that describes, in a probabilistic fashion, a physical event. The procedure by which this is done follows from the fundamental principle of Monte Carlo (16 ), which may be stated as follows. If p(x)dx is the probability of x lying between x and x+dx, with a < x < b -b p(x)dx = 1 (41) a Then r = P(x) = j p(y)dy (42) Ja determines x uniquely as a function of the random number r which is uniformly distributed on range (0, 1). The quan- tities p(x), p(x) are the probability density and the probability distribution functions, respectively. With the development of large computers, the use of the Monte Carlo method in a wide range of problems has 37 increased rapidly. In particular, the method is applied to various physical problems such as plasmas (15), radiative transfer problems (21), neutron transport and nuclear reacter physics (29), etc. The random walk type simulation can be briefly described in the following manner. One first introduces the basic physics of the problem into the computer in a probabilistic fashion. A system of coordinates and boundaries are defined and then, as a computer experiment, particles are released from the source. These particles are traced as they diffuse through a prescribed medium following the probabilistic interaction laws. The particles are followed until they escape from the medium or are absorbed or are converted to the thermal field. The parameters pertinent to the evaluation of the desired quantities are recorded. One continues processing additional particles until adequate statistical estimates of the quantities of interest have been obtained. When doing a Monte Carlo problem one focusses attention on three main topics. They are: 1. The development of the analogy for the probability processes (simulation). 2. The generation of sample values of the random variables on a given computing machine. 3. The design and use of variance reduction techniques. 38 The variance reduction methods are often strongly dependent on the probabilistic model, thus the greatest gains are often made by exploiting specific details of the problem, rather than by routine application of general principles. B. Importance Sampling The general idea of importance sampling is to draw samples from a distribution other than the one suggested by the problem. Then an appropriate weighting factor is introduced to correct the biasing caused by changing the original distribution. If done correctly the final results are essentially unbiased. The object is to concentrate the distribution of sample points in the parts of the region that are most "important" instead of spreading them out evenly. Thus, the probability of sampling from an interest- ing region is increased; the probability of sampling from unimportant or less interesting region is correspondingly decreased. Stating this formally, we have e = j% 1, V- then el i ' ' 1 2|' E> 24l (E ) i 77^ ~T£ro "i = Kl (53) and c « w i- l.. C i-TT- (i ^> (53a) 1 el Equation 53 holds for low energy region, which is approxi- mately independent of energy. • C. Russian Roulette and Splitting In general, the sampling is done such that we examine the sample and classify it as being in some sense "inter- esting" or "uninteresting" . We are willing to spend more than average amount of work on the interesting samples and on the contrary, we want to spend less effort on the uninteresting ones. This can be done by splitting the interesting samples into independent branches, thus resulting in more of them, and by killing off some of the uninteresting ones. The first process is splitting and the second Russian Roulette. 43 The "killing off" is done by a supplementary game of chance. If this game is lost, the sample is killed; if it is won, the sample is counted with an extra weight to make up for the fact that some other samples have been "killed". The game has a certain similarity to the Russian game of chance played with revolvers - whence the name. For the present plasma system we divide the whole energy spectrum into two regions, the low energy region I (0 eV to EDI eV) and the high energy region II (EDI to E ). Here E_ is energy of source particles, which is the highest energy the particle can have. EDI is the energy region dividing line, which can be chosen as any value. between 5 eV and 100 eV depending on the.initial guessed energy distribution. In region I, the number density is high, and the elastic process is the dominant process. As we mentioned earlier, the elastic process does not contribute much, so Region I is classified as uninteresting. Each particle in this region (particles with energy less than EDI) is assigned a weight kl(kl > 1). In other words, each particle plays the role of kl particles. . In region II, classified as interesting region, each particle carries a weight w ? (w ? = /k ? , w ? < 1). Namely, each particle in this region is split into k ? particles with weight w~. As particles cross over to other regions, care must be taken to account for killing off and splitting effects. 44 For instance, a particle in region I may acquire enough energy, either from electric field or by superelastic collision with the excited atoms, to cross over to region II. Then we should generate k~ particles (k~ = k. • k ? ) in region II with the same energy. Particles from the high-energy region may lose part of their energy by inelastic collision such that they cross over to the low-energy region. Then killing off should be effective. For every k~ such particles that cross over from region II to region I, only one particle survives. The values of k, and k~ can be chosen according to the total number of particles in the system to be processed, the source rate, and related accuracy requirements. They are input parameters. D. Initial Guess Monte Carlo methods comprise that branch of experimental mathematics which is concerned with experiments on random numbers. The user, like the experimental physicist needs theory and knowledge to give structure and purpose to his experiments, and as experimental work provides growing insight into the nature of a problem and suggest appropriate theory. Good Monte Carlo practice should keep this relation- ship as a general maxim. The basic procedure of the Monte Carlo method is the manipulation of random numbers, but these numbers should not 45 be employed prodigally. Each random number is a potential source of added uncertainty in the final result. Thus it will usually pay to scrutinize each part of a Monte Carlo experiment to see whether that part can not be replaced by an exact theoretical analysis contributing no uncertainty. In other words, exact analysis should replace Monte Carlo sampling wherever possible. Sometimes reliable intuition would aid in increasing the efficiency of Monte Carlo calculations. For the plasma system we are simulating, the straight-forward approach would be to follow N source particles on their life time history, until they lose energy sufficient by collisions to reach a thermal equilibrium. Over 50 inelastic collisions generally occur, and the probability (inelastic cross sections) for inelastic processes is low. Consequently, the particles travel a relatively large distance which consumes a good amount of time. Thus a large number of time cycles would be required before the system reaches steady state. To avoid this, we use our intuition, to start off with a guessed initial distribution. The electron energy distribution f(E) (or electron number density distribution) is a function of the electron's energy E. Initial particle data is generated using this distribution, and this can save considerable irrelevant computations. This idea is simple, yet very useful and effective. Without an initial distribution, the problem would require hundreds of time cycles to 46 converge. However, with a reasonable guess, good results have been obtained in as little as tens of time cycles. Of course, the better the initial guess, the faster the convergence. E. Antithetic Variates When we estimate an unknown parameter 9 by means of an estimator t, we may seek another estimator t 1 having the same (unknown) expectation as t but a strong negative correlation with t. Then l/2(t + t') will be an unbiased estimator of 8 , and its variance is Var[-|(t+t' )]= ~var(t) + jvar(t') + icov(t,t«) (54) where cov(t,t») is negative, and var(t) = £ [(t-u) 2 J cov(t,f) = £ f (t-uMt'-u«)) u,u' are the mean of t,t' respectively, and£(x) is the expectation of the random variable x. It is possible to make var'(tft' ) smaller than var(t) by suitably selection t'. For example, 1-r is uniformly distributed whenever r is, and if f(r) is the unbiased estimator of 6, so is f(l-f). When f is a monotone function, f(r) and f(l-r) will be negatively correlated. Thus we could take ■|(t+t« ) = \ t(r) + -| f (1-r) (55) 47 as an estimator of 9 . The name and idea of antithetic variates was first introduced by Hammer sley and Morton (13, 14) during attempts to improve Buffon's needle experiment for estimating the value of n (3.1415926) . The main idea of antithetic variates is based on the underlining theorem [from Hammersley and Morton (12)] which can be stated informally as follows: "Whenever we have an estimator consisting of a sum of random variables, it is possible to arrange for a strict functional dependence between them, such that the estimator remains unbiased, while its variance comes arbitrarily close to the smallest that can be attained with these variables". As the name implies, antithetic variates are the set of estimators which mutually compensate each other's variations. Essentially, we rearrange the random variables by permuting finite subinter- vals in order to make the sum of the rearranged functions as nearly constant as possible. From the practical viewpoint, the mathematical conditions imposed by antithetic variates to the Monte Carlo calculation are quite loose and flexible. It is relatively easy to find a negatively correlated unbiased estimator to produce an efficient variance reduc- tion scheme. A system of antithetic variates is obtained by simple transformations based on a stratification of the interval (12). The transformations are constructed in such a way 48 that the efficiency of the method will increase as a higher power of M (the number of antithetic variates taken). The transformations considered are linear combinations of the values of the function f at a number of points. Thus a correlated stratified sampling technique is developed, and its efficiency generally exceeds that of crude sampling by 3 a factor M . The following useful transformations are obtained in this manner: t. = - T\ -faf(ar.) + (1-a) f[a+(l-a)r.]} 1 n *— J , <• l l * i=l (56) n n <-■ n a i i=l 1 n 2 = n f± 1 {af(ar. ) + (1-a) f [l-Q-a)^]} (57 t„ = n l-l The value of a is between and 1, the variance var(t.) has a minimum at some value of a. An adequate rule of thumb is to choose a by finding the root of f(a) - (1-a) f(l) + af(0) (58) Another useful transformation is ^ m t _ = if f(£±l) = U f(r) 3 m 4-" ^ m m m r=o (59) 49 In general, the variance of crude Monte Carlo is var(t) = 0(m~ ) while the variance of t~ is var(t-) = — 2M 0(m ), (M "> 1). The simplest form of antithetic transformations have proved quite successful in neutron transport calculations. The general notion and a simple form of antithetic variates are applied in our present calculations. High efficiency gain can be obtained easily if it is used sparingly and at judicious points in the computation. We have used the random number r.(r. e(0,l)) most of the time for playing the law of chances, and 1-r . is another random number of the same kind which is the image of r.. The two are used together in good many places where repeated use of r. is required. This serves the purpose of two antithetic variates and cuts the labor in half (same random number r. being used twice). Equations 55 and 59 have been used in the initial data preparation program and source routine where initial particles states (E . ,u . , z . , wt . ) , E. are generated according to input initial distribution. For instance: u. = 2r.-l, u i+1 = 2(l-r.)-X z i " r i d ' z i+l = '<» E i " \ + r i- iE k' E i + 1 " E k + -A E k <60) E . = E n for source particles. Here k is an energy index, or kth energy interval. Wherever r. appears, r. and 1-r. are used interchangeably. 50 Thus we are generating source particles from diametrically opposite directions, or from symmetric locations. F. Double History (Double Processing) Technique In the particle tracing routine, using the general idea of antithetic variate technique, we have adopted a "so-called" double history processing technique which is described as follows; a) In the course of tracing the particle's history, in each time cycle, we use the transformed random number X* I lnr to decide the distance to collision d n , and go 1 col 1 r on until the particle reaches census time. However, we process the same particle twice under the identical initial condition except that we use >^«|ln(l-r)| as the distance to collision d ., „. Both sets of history are col 2 J ta lied but only one set of particle data is 'stored. This way each particle is processed twice independently but correlated in such a way that they compensate each other's variation (variance). b) Another alternative for double history processing is, to play a dual direction random walk instead of dual collision distances. We process the same particle twice in such a way it looks like there are two identical particles emerging into the system in two diametrically opposite directions. Again both histories are tallied, but only the original particle's data is stored. 51 The two versions of the double history processing technique increase the labor (computation time) by a factor of 1.2 to 1.5. The efficiency gain is difficult to specify, At least it is certain that the gain is more than double the total number of the particles in the system. The purpose of dual processing is to compensate each particle's variation and minimize the variance (statistical errors). One set of particle data in dual processing is thrown away, while the other set is stored just like normal processing. Which one of the two resulting sets is saved is completely random. 52 CHAPTER IV MONTE CARLO SIMULATION CODE A. Introduction This chapter piesents the complete Monte Carlo simula- tion program, including an input, output, main flow chart, and the background theory and techniques incorporated. This program implements the condensed history approach that was described in Chapter II together with various variance reduction techniques described in Chapter III. The basic Monte Carlo calculation can be described most easily by reference to the simplified flow chart shown in Figure 4. By means of this flow chart we shall follow the particle tracing process, making relatively brief comments on the various steps implied by each box, on each routine, and on the main program variables. The program is written in FORTRAN IV and is designed to be run on IBM 360/75 under the HASP-MVT system. B • Calculational Procedure and Program Description To obtain an overall view of the calculations, we note that there are three basic steps involved in the simulation process. The first is the initialization of all parameter values, i.e., prepare the particle data and predict the unknown dominate parameter values. The second is the main loop of particle tracing which generates census information. The third uses the census information and takes tallies to estimate the histograms (electron energy or number density 53 Input control ' * Initialization and particle data generation (E, Z, U, WT) < ' Follow initial source particles to census, tally particle states i " last source particle? No * Yes Follow residue particles to census and tally states i . 1 f last residue particle? No • Yes From state tallies (particle histories) calculate estimatd parameter values + Renormalize and correct the ion density, reinitialize parameter values increment t n to tn+1 time pycles and print the output quantities > ' If a prescribed number of time cycles have been processed, output and stop. Figure 4. Flow chart of the Monte Carlo simulation code. 54 frequency distribution vs. energy), and parameter values. Finally the dominate parameter va ues are corrected and the program returns to the first step. Let us follow the logic flow diagram (Figure 4) in the order of the number marked on the left hand side of the boxes, and describe their individual functions. 1. Input Control: The main program is usually compiled and the loaded version is then stored on disk storage, ready for execution. In order to offer a wide flexibility, the following input parameters are used to control the program runs. a) NDISK : program start control. The program can start from the beginning (cycle 0), or it can start from the stopping point of the previous run. This feature enables long runs to be broken up into several short runs. It also offers debugging flexibility and economy of computer time. While a long run may run up to 150 time cycles, we can always break it up into multiple numbers of fives and tens cycles. Thus we can adjust other parameter values on a cut and try basis without wasting much computing time. The input parameter for this control purpose is NDISK: input for storing from very beginning; input 1 for starting from the latest tenth cycle of last 55 run; input 2 for starting from the latest fifth cycle. b) Kl, K2, and EDI: Regional weights and dividing energy. The Russian Roulette and splitting variance reduction techniques use three input parameters to control the amount of killing off and/or splitting of particles. This also determines where to divide energy regions for the Russian Roulette and splitting, assuming that the option to use Roulette is selected. The input parameters are Kl: the weight of particles in region I (Kl '■ 1) , K2: the inverse of the weight w2 of particles in region II where w2 = 1/K2 (K 2 '• 1) , EDI: the energy dividing line, namely (0 to EDI) is region I, (EDI to E ) region II where the energies are in electron volts. c) NTC : The total number of time cycles to be run. DTC : The input parameter to control the length of census time At of the time cycles. In terms of distance travelled by the particles, At is equal to a fraction or multiple of one mean free path , corresponding to average energy of particles in the system. 56 d) NTP: The total number of particles that actually exist in the simulation system. NTW: The artificial* total number of particles in the system during Russian Roulette and splitting games. NS : The number of source particles injected into the system during each time cycle. e) KELC: The input parameter that controls the impor- tance sampling option of the program, i.e., the multiplying factor used to suppress the elastic cross section (reduces the probability of elastic process by a factor KELC). f) EDP : The input value of the ratio of the electric field strength to the pressure. PDPO: The input value for the reduced pressure of the system ( P n where p« = 1 torr). g) NIP: The input parameter used to control double processing: NIP = 2 for double processing, NIP = 1 for single processing. IUSN: IUSN s 1 for dual-collision type double processing IUSN = -1 for dual-direction type double processing *The total number modified by scaling and weighting factors. 57 h) D: The input value of the width of the slab (diameter of the cylinder) which defines the external boundaries of the medium. JZ: The number of zones the slab is divided into, or the number of internal zones. i) NCUM: Inp t control to give cumulative results (time average) over the time cycles or a cumulative average over NCUM time cycles. NAVP : Input constant to control the power of At, (or &D) namely At.**NAVP is the weight of the cycle to be used for cumulative results. These two parameters are used to implement antithetic variates techniques along the time axis. j) EO: The energy of source particles. SR: The rate of the source particles entering the system. k) FID: An array that contains the initial guess distribution (number density vs. energy). 1) NUSC: An input parameter used to control the up- scattering elastic collision process; NUSC = 1: with upsacttering, NUSC = 0: without upscattering. 2. Data preparation: The initialization is entered only once during the calculation, where the time cycle is (or at time t„). After the input phase, all the pertinent 58 parameters are set to their proper permanent or temporary values. Before entering the next phase to trace the particle histories, several dominate parameter values must be predicted. The data preparation routine is called in to generate particle data (E.: energy, u..: direction of travel, Z.: position, Wt . : weight) based on the guess for the initial distribution (generate machine particles). Then the average energy of the system is calculated and is used to predict the census time (At or equivalently DCEN) as well as the number of ions(N + ) . 3. Particle tracing: With the point of origin, direction of propagation, and energy (velocity) known for the initial source particles, the next phase follows the particles through a series of stochastic interactions until the census time is reached, and the particle data and parameters values are tallied and updated. During the course of tracing the particle histories, all the relevant events are recorded and later on can be inter- rogated to provide the output quantities of interest. This section of the program is the central part of the main program. It contains the main loop for tracing all the particles in the system in a self consistent way. Most of the variance reduction techniques are incorporated here, including double processing, Russian Roulette, and splitting, as well as the various versions of 59 importance sampling techniques. 4. Secondary particle tracing: The residue particles, i.e., particles generated by ionization collision process (the secondary electrons so produced) are treated in all respects like the regular source particles, except that their starting positions and time coordinates as well as directions are those associated with the particles that generated them. These residue particles may go further and produce other residue particles (secondary electrons) if their energy is high enough, thus forming a short production chain. The probability of such chain produc- tion is small however, since the average energy of the secondaries is quite low. The end of the main loop comes after all such residue particles are processed. 5. Data calculations: The details of the phase of the computations depend upon the specific goals and charac- ter of study. Our primary interest is the particle energy distribution (or number density vs. energy). This frequency distribution and other pertinent para- meters are estimated (or calculated) using the stored data from various parameter counters in the main loop (tallied information) as well as tallied histories. The quantities of interest include local values (for the current time cycle only) and the cumulated values from several time cycles. After these calculations, the next phase is to prepare for the processing of new time 60 cycle. 6. Normalization and printout: Up to this stage, all the output quantities have been calculated. Next the output information is printed out, and at the same time, if the cycle number is the multiple of five, all the particle data and related parameter values are stored on the disk. Then the next run can be read off of the disk and the process can proceed from where we left off the last time. All the parameter values are again initialized, and new source particles are generated for the next time cycle. Particles that have escaped from the system or have been absorbed create holes in the particle data storage which are filled up by the new source particles and the residue particles produced during the time cycle. Before going on to the next cycle, the ion density is renormalized and corrected in a semiadaptive way accord- ing to a balance factor BF. The ion density distribution parameter N is a dominate parameter since it affects several processes such as the rate of recombination and the rate of super-elastic collision (refer to Equations 27 and 28 ) . It is corrected in such a way as to accelerate convergence, till the system reaches a bal- anced situation. Then increment time from t to t . The output quantities are the following: a) The kind of gas in the medium, source particle energy, pressure, slab thickness, and At the time interval of this cycle. 61 b) The number of time cycles, the time in seconds, the number of source particles generated during the cycle time, and the average energy of the particles in the system. c) The number density distribution vs. energy, the single cycle distribution, the cumulated distribu- tion, the number-density deviation, the average- energy deviation, the average energy, and the energy interval vs. energy. d) The number density and the average energy vs. z (position zone) . e) The ion distribution (N+), the number of ionizations (NIZN), the number of escaped particles (NES), the number of absorbed or recombined particles (NAS), the number of super-elastic collisions (NSE), the number of excitations to metastable levels (NMET), and the number of excitations to other levels (NEXC). f) All the input control parameters as listed in section IV-1. g) The electric field strength (EF), the source rate, the W-value, the average energy of escaped particles, and the average energy of absorbed particles. h) The number check on the density disbributions. All the above outputs include both current cycle values and the cumulated values over the past cycles. 62 Final checks: The last stage is to check whether the prescribed number of time cycles have been processed, adequate statistics have been collected, and all the useful data have been stored on the disk. If so, the program stops. Further programming details are provided in Appendix 63 CHAPTER V RESULTS A. Convergence and Reliability The present Monte Carlo Model is a predict-correct iterative scheme, with semiadaptive control ability. For such an iterative procedure , our main criteria for the evaluation and determination of convergence are based on the following inferences: 1. The balance factor (BF) as defined in Chapter II (Section II-B-3) signifies whether or not the system has reached steady state, i.e., converged to a final solution. At steady state, the value of BF should be approximately equal to unity. It may approach to unity in both directions, and due to semiadaptive control procedure used, it may oscillate around unity. Larger amplitude oscillations are observed at first, but they rapidly damp out. 2. The estimated values of the major parameters, e.g. the average energy and the W-value show a similar fluctua- tion, and then they converge to nearly constant values as steady state is approached. 3. The validity of the convergence has been checked using test problems where analytic results are available. A good comparison has been obtained, showing that the program and the model are working correctly, free of obscure errors. 64 B. Error Estimation and Analysis 1. Error Analysis 1. As mentioned in early chapters, Monte Carlo simula- tion is equivalent to performing mathematical experiments, the results being based on observations of such experiments. From this analogy we can derive the following procedures for the analysis of errors or accuracies of our calculations. It is convenient to subdivide experimental errors into three broad types, namely: random errors, systematic errors, and blunders. In general, the experimental error is some additive function of all three, while blunders can hopefully be eliminated. We shall describe them in detail in follow- ing sections. a. Random errors - such errors are of great concern in Monte Carlo type calculations. They generally involve statistical fluctuations or deviations, representing the difference between the singly measured value and the best value of a set of measurements, i.e., the difference between the arithmetic mean and the true mean. Sometimes such errors are referred to as ele- mentary (or inherent) errors in the measurements. The commonly used measure of random errors is the variance or standard deviation. b. Systematic error - a systematic error tends to have the same algebraic sign, ± . e. it' is either an additive or subtractive quantity introduced in the measurement 65 process. It is an unpleasant and insidious contribu- tion which is not generally amenable to statistical treatment. Thus such errors seriously impair the reliability of the estimation. Typical examples include: i) Incorrect assumptions or approximation in the representation of certain processes, ii) Constructional faults and mistakes in the algorithms or subroutines, iii) Inadequate regard of constancy of experimental conditions or inadequate sampling techniques (biased) . c. Blunders - these are outright mistakes which should be corrected by all means. Possible examples include: i) Incorrect logic, misunderstanding of the problems, ii) Errors in transcription of data, iii) Mistakes in constants used, iv) Confusion of units. In any type of calculation, the systematic error as well as blunders should be removed. One way for detecting such errors is to compare the program against some known reliable solutions. Only the random errors are subject to reduction by the various treatments amenable to attack by variance reduction technqiues. The precision in an estimated (mean) value is propor- tional to the reciprocal of the statistical error and is 66 high if the statistical error is small, i.e. the accuracy is high if the net systematic error is small. Usually, but not necessarily, high accuracy implies small statistical errors as well. As discussed in later sections, our model and program has been checked against several test problems, and good agreement has been obtained. Therefore, it can be assumed that any serious systematic errors have been removed and that the program is free of blunders. Thus, we shall con- centrate on the analysis of random errors. No thorough analysis of variance has been made. Due to the intricate way that the histories calculated during a given time cycle depend on previous time cycles, a meaningful variance estimate beyond the first time cycle seems to be out of the question. However as the system reaches steady state, the deviations or variations (statis- tical fluctuations) of the estimated parameters can be estimated by means of standard deviation (or variance) over subsequent time cycles. Namely, we use the following formulation : S. = 1 £. (f . . - f .) n V i (61) where s . : l f . ij f . : .i the standard deviation at the ith energy interval over n consecutive time cycles. the parameter value (density distribution func- tion of ith energy interval, at jth time cycle. the average parameter value of ith interval over n time cycles. 67 A heuristic approach was adopted to measure local deviations or statistical fluctuations at each time cycle. This pro- vides a relative magnitude of deviation or dispersion which can be used to compare the various variance reduction tech- niques. The approach is described as follows: The estimated major frequency distribution (the elec- tron energy density distribution function in our case) is obtained in the form of histogram. The abscissa of the histogram is divided into so-called class intervals (energy intervals in our calculations). In each such interval an erect rectangle (block) of heigh f. and width E. is formed. 1 1 This block-area type of distribution is the fitted fre- quency distribution curve. We devise a numerical descrip- tion of local deviations by defining: i) the location index of the center of the erected rectangles in the histogram as i and with abscissa value E . . l ii) the spread or dispersion around the center. The total number of particles in each class interval (energy interval) is proportional to the area of the rectangular flock (f. AE.) / the particle's energy in this class interval falls in the range (E. - AE . ._ f E. + AE. /? ) with overall average energy E~. = E.. Then the following relationship holds (refer to Figure 5): 68 - k UJ AE, «*- i > AEj*| o c Z3 , . c o »l f M 3 5 < ' i ' ^_ E| E i + I E l+2 E-Energy, eV Figure 5. Rectangular block representation of the electron energy distribution histogram. 69 ifi til AE. " f . AE or A£ i - f i AE7 (62) where AE~. is the deviation of the average energy FT. from E. or AE~. = E\ - E.. Also Af. is the corresponding deviation 111 1 r -a of the frequency distribution at location index i. Equation (62) provides a relative measure of the ampli- tude of the deviation of density distribution curve at energy E. which can be used for comparing the efficiency of different variance reduction techniques as well as provide an heuristic clue to the modification of the fitted density distribution curve. 2. Data Filtering The random errors (statistical fluctuations) result in bumps and ridges in the estimated frequency distribution function. It may be possible to alleviate many of these distortions by means of data smoothing techniques (data filtering) . Two types of data smoothing techniques commonly used in the processing of time series might be used here: a. Moving averages - this technique operates by replacing each point of the frequency distribution function (height of the histogram at E.V.) at time cycle t. with an average value of several subsequent points in the time cycle series. Thus if f . . is the value 70 of frequency function at ith energy interval and jth time cycle, then f . . is replaced by J/M where W., the weighting factor for the jth time cycle, usually has the value of 1 or we can define it as a function of &t . . The value of M is the number of successive points of the time cycles to be included in each average, b. Parabolic filtering - this is similar to the moving averages technique, but M is selected as an odd number, the point to be replaced is located in the center of these M points, and a parabola is fitted to these points by means of least squares. Thus it is replaced by the corresponding point f . on the derived parabola. The process is repeated by shifting the center point to the next one in the time cycle series. The technique of moving averages (a) is employed in the program, because of it's simplicity and the fact that the parabolic smoothing technique is too sensitive to the changes in the distribution function. Note that this smoothing (filtering) technique only becomes effective as the simulated system reaches steady state. 71 C. Solutions of the Problem 1. Preliminary Checks ■ In order to verify that the model is valid and correctly working, and the computer program is free of bugs f blunders, and systematic errors, we solved a specific problem first reported by J. A. Smit (11) and later extended by Heylen and Lewis (10). Their solutions for the electron energy distri- bution were obtained analytically based on the Boltzmann transport equation for the special case characterized as follows: i) electric field present ii) no external sources iii) infinite size medium of background noble gases. The comparisons are made for helium gas at E . =4, 5, 10 V-. cm~ . (E^_ . : electric field strength to pressure ratio), ' f/p The results are shown in Figures 6, 7 and 8. Both Smit and Heylen' s results have taken into account all the collision processes including those due to elastic loss as well as inelastic losses, Heylen' s results is more up-to-date, for lowE^., values onlyE,. ', = 5 is given, which was chosen for f/p J f/p the present comparison. In so far as possible, the same cross sections as Smit and Heylen and Lewis used were incorporated in the present calculations. The slight discrepancies apparant at the low-energy end is probably due to differences in treating the scattering process, and 72 u o tt V CD IH CM w ■p X •H o> M >% •H • <0 0> T3 VO c l±J m|w in ' Q) ^ • UJ Q) ^ CJ VO II lectron on (193 ion (1? GO of the e alculati calculat ° e (3) J 'uojiounj uojinqujsiQ •H lO 73 CM fO O fO 00 in II 10 CM W ■p CM CM M CM O m O o CJ > •H 0) -P *-» 3 H CO ,Q in •H 0> -P<- c W K M O 4-1 C O •H *J XI •H M ■P (0 •H TJ II M 0) .— C O II n, •P w O w •H 4J (DHU 3 H m w nj U C W •H U id w CO 0> •H (3)1 'uojpunj uoijnqujciQ 75 the many approximations required in their analytic solutions. Also, some of the cross sections used in this calculations are experimental values, which may differ by some amount from that used for analytic solutions. For the high energy region, (near the source energy) checks have been made against the analytic solutions reported by Lo and Miley (32). As shown in Figure 9, the general shape of the electron energy flux distributions are in reasonable agreement. In this case it was not practical to use the same cross sections as incorporated in the analytic solutions, and this may account for some of the differences observed. Also the analytic solutions are obtained for an infinite medium vs. a 1 cm slab in the present case. Thus the analytic solutions do not allow leakage, although this becomes important in the present calculation for low pressures. Lo and Miley' s results are plotted in Figure 9 for the source energy E n = 1500 eV and E = 500 eV, whereas the present results are for E n = 1000 eV. The plot shown is for the flux density (normalized) instead of energy number density. They are related as follows : CyP(E) = f(E)-V(E) (63) where ^(E) f(E) V(E) electron flux density electron number density Electron kinetic velocity 76 _t-t t r 1 'I — T -T- MM — 1 1 1 1 — - > 4) 8 ii o C — — > O) O o ■"* II \. — _ UJ° > 0) O O ifiC ii 3^^ <] \ ~— Ui / T - M J Ml 1 1 1 III I 1 1 1 1 O O > M 3 -P O H P >i W +J CD •H M (0 C W i a> X H 3 -H .H S m "> T) a> I ^ 1 >% C ( o> ^ •r i_ •H 4 a> ■p c UJ > M in r 1 1 CD i C II 1 0) 1 1 o 1 ,C W 1 0> : ■H«0 X! C fO &> 4J 0) m o ( o o ( in c C H i o W II •H H O f0 w I ft e o~ ' U (0 i o CD co o> v. c Ixl Id "2 (3)^uo|punj uojinquisja a> •H 4-1 O TJ m a) ■P O > O 0) H o o o c o nJ h ■p II P o jC W ■p •H » ? M M 0) K P ^ O H m II C 1 o &< •H -r P -P 4 3 ro : X) ^ •H W •- M d) ) -P P -1 w rt i •H )-i •. T3 n >i C P 3 r +Jr c c o » +j c w 3 Q> o d) ,Q c •H -P -H UJ •P 03 M P ^ -P 9-. ^-» .Q W o •H Q) -H UJ rgy distr ent sourc lus 1/E d T-L 0) M ft o electron e Id but diff Maxwell i an M (U (D-O o Xi -H c En m 03 (3H u 0!P un d uoiinquisiQ Q) •H o^ o «4 o - © U) C LJ •H * *— o o > E 0) o c r* >s ii • o> c c 1~ o a a> •H -H c •P +> +J bl 3 «I 3 ,Q A °c> ■H W -H LJ ergy distr ource rate 1/E distr — c w w b 0) 3 electron different wellian pi M 0) 4J X 'o 43 3 «J E-< ,Q S (3) J UOipunj i{oj;nqu|S!Q m H o c UJ T3 a) •H m o •H M +J o .Q •H M ■M W •H TJ C rtJ •H H H Q) X rtJ s M C fd -P w a) ■p x: •p •H u a> w I ,1 2 C m O E W •H \ M # fO o o U H (3)j uojp'snj uoj|r,qu|sjQ a) u •H PL. 85 -O o o o •H <+H O •H M •P O K W h o Oi C H v. o a> •H II c 4J o UJ 3 o •— * X) W UJ ergy distri ressure at o c a electron e different 'o (D H > P m (3)) uojpunj uojinqi j|sja in U •H T3 (3) J. uoipunj uoj|nqjj;sia H CD >» a> V- Q> C UJ 'o • O 'O o 2 cu c o •H P 3 X! •H M -P co •H •0 W P M o p o VD p II CD c o P p o CD ft cd m w CD £ p B CO (3)JU0!pun_j uojinquiSjQ CO H CD u •H 94 o o O o g «> g g II H n ii o CO o CO o CO o CO I _J L 2 M O °o > >s a> c iLl rH 0) •H — c> 0) •H •P o P r* .0 ■H II M •P w w •H T3 •» O >iH CP M II K 9 M « C o > "~ — m a >, o» c o la O r- « •H c -P II u 3 ,Q o h ""■^ •H W UJ nergy distr 760 torr, o e electron e /p - io. P - '2 H W « O CN 0) M 6 •H 96 1 i i 1 \ \ 1 / / / — r i \ 1 f 1 1 / 1 J 1 1 1 \ 1 \ \ • ■* \ — %> Q Q \ H it \ \ o o \ \ v> en \ \ i — \ \ 1 \ \ i 1 \ \ 1 \ \ . \ \ \ \ ! l\ 1 \ I ! 1 1 I | - '2 u u o •p o & - % w o o \> •o b b 0* \ m W ■P «J c • <— » •H *-» > ■P 0) C W c p 1 m n E >» c u 9 \ w i W 0> >H TJ a) c c *0 rH ffi © O C m >» c o> k_ •H O O c ■Q •H LU M ■P 4 (0 > •H CD TJ O & O O M H (1) 1 C II O c ( u •» +J kO u H 0) O H H

'O w (3)J uoijounj uojinquieiQ CM CM CD U & •H P4 98 10 O o — 2 ^-. UJ (0 •c H 0) •H HH U •H ■P U 0> C •H 5 0) > u o m c w >» c o •H ■P tf CJ r 0) H 0) 0) g ■ m (3) J uoijounj uojinquisjQ cm 9 M 3 CP •H h 99 > ill UJ +J id 4) u o m c o •H +» ,Q •H M •P • W > •H Q) T> O >iO D^O M H •P M U M O QJ N || •H rH ft (0 e - M O O C 0) II ft m (3)_N uojjounj uoi|nqu|Sja CM a) •H Cm > 0) ( C I P -i •H •> s : P 'rd ; (3)) uoijounj uo!inqu|siQ in £ O En »W OH uoipunj uojinquis.iQ M •H P4 106 2. The dependence of the distribution functions (with or without electric field) upon the pressure and source rate only shows up in the high-energy region, where the distri- bution functions have higher tails for higher source rate, and have smaller values for higher pressures. 3. The range of variations observed in the W-value . calculation reflects statistical fluctuations which limit the calculational accuracy of the absolute numerical W- values. Still the results demonstrate that the W-value is more or less independent or insensitive to the changes of pressure, source rate and source energy. 4. Due to a tight budget, the primary emphasis was on increasing the computation efficiency and reducing the computation time to a possible minimum. Thus most of the calculations were carried out with 2000 to 5000 actual particles and for any particular set of parameter values, and the results were generally obtained within 10 minutes of (IBM 360/75) computer time. Compared with hours of computer time and 5 to 10 times more particles used by others for the similar calculations (e.g., see References 3, 4, and 9) the gain in computation efficiency seems to be significant. Of course, a complete comparison is not possible since all of these calculations achieved different accuracies in the final results. 5. The present results also show close agreement with the analytic solutions of Lo and Miley (32,33) on the 107 following points: a. The shape of high-energy portion of the distribu- tion curves are in reasonable agreement as shown in Figure 9. b. The high-energy part of the distribution function is not changed by the presence of an electric field for field values E_ , < 10 . f/p - c. With an electric field, the low-energy part of the distribution function agrees with Smit's result; and this part is not affected appreciably TO T by source rates < 10 particles/cm -sec. 108 CHAPTER VI CONCLUSIONS A. Prelude The application of Monte Carlo techniques to solve particle transport problems is of fundamental interest in many fields of physics and engineering; especially when analytic and/or numerical solutions of the basic governing equations are too complex to be practical. Yet direct application of Monte Carlo techniques to simulate individual interactions is sometimes prohibitively costly because of the large amount of computational time required. Thus techniques for improving the efficiency are required. Various existing variance reduction techniques can be used, but they are problem dependent. Thus techniques that are efficient for one type of problem might not be as effective for other types of problems. Thus such techniques must be judiciously selected and modified for the particular problem under consideration. One basic underlying principle which applies to any Monte Carlo calculation is: Apply as much known information as one can (given in analytic and/or numerical form) to reduce the uncertainties of the problem. Whenever possible Monte Carlo experiments should be checked and replaced by exact theoretical analysis to reduce uncer- tainties . The Monte Carlo experimentalist has to exercise ingen- uity in distorting and modifying problems in the pursuit 109 of variance-reduction techniques. Although Monte Carlo methods are general devices, still a great deal of work depends upon the individual's originality to create special methods or numerical schemes to suit their needs. Unlike the physical experimentalist, the mathematical experimentalist using Monte Carlo Methods has the advantage that his experimental material consists of mathematical objects which can be distorted, controlled, and modified more easily. The greatest successes of the Monte Carlo method have arisen where the basic mathematical problem itself consists of the investigation of some random processes. However, there are exceptions involving deterministic problem- solution of boundry value problems and partial differential equations. The solutions of these problems are closely connected with the characteristics of certain random diffusion type processes (or can be converted to such type of processes), therefore these problems are reduced to the modeling of such processes. Thus the model and the special techniques developed in present research are quite general in scope, and it should be possible to apply them to many classes of Monte Carlo calculations with a certain amount of modifica- tion. B. Summary of Present Work In this research study, we have devised a general mathematical model for the particle transport or diffusion 110 type plasma problems which have inherent nonlinear proper- ties due to recombination. Many special variance reduction techniques are incorporated into the system. The model can be used for calculating both time-dependent and steady-state nonlinear transport problems provided that suitable methods exist to calculate, predict, and correct the important parameters which dominate the nonlinearities. A wide range of nonlinear problems can be formulated and handled in this way. The principle idea, the piecewise linearized predict- correct model, is a general, simple, and efficient approach for solving nonlinear problems. In principle, most nonlinear problems can be cast in terms of the basic ideas and algo- rithms of present work. This represents an extension of earlier works of Musha and Itoh (8), Thomas (9), and Fleck Jr. (11). As for the amount of computation time and effi- ciency involved, the present study has made significant improvements. Some of the improvements are a natural con- sequence of the large storage space available in modern computer systems. Additional improvements come from a few simple, but effective, special techniques incorporated into the present model. As described in Chapter II, these include: 1) The initial distribution - straightforward particle tracing procedures would require following every simulated particle from its birth (from source) with initial energy E at t = t to the end of its life time 3J o o history. This type of simulation would require hundreds Ill of cycle times to achieve a final steady state condition. Instead, we start with an apprixmately steady-state system by guessing an initial particle energy distribu- tion. Then we simulate the system by tracing particles starting with energy corresponding to this input distribution. This is eguivalent to skipping over the initial transient time cycles and greatly reduces the computation time and increases the calculation efficiency 2) The application of the negatively-correlated variates (antithetic variates) technique greatly reduces statistical fluctuations with a minimum number of simulating particles (machine particles) . 3) The introduction of a weighting scheme and an artificial game of chance (using weighted cross sections) contri- buted some reduction in calculation time and resolved some practical difficulties. 4) For the system with an applied electric frield, the exact formulation for a parabolic flight path was used to advance the particles in larger steps than would have been possible using the small-step linear approximation employed by Musha and Itoh (8) and also by Thomas (9). This improves the efficiency appreciably. Throughout the entire calculation, one basic balance equation (Eq. 30) played a very important role. It merely states that at steady state, the source is balanced by losses 112 By means of this equation, we predict the nonlinear parameter values and provide checks for convergence of the solution. The need for large amounts of storage space for Monte Carlo calculations no longer creates serious problems or drawbacks on modern large computer systems , such as was the case as recent as seven or ten years ago. In addition, techniques exist for further reducing the storage require- ment. In the present study, we suggest a technique to eliminate major storage requirements, and use only a small amount of storage as a temporary buffer storage (described later) . In other types of numerical solutions, e.g., in finite difference solutions of the diffusion equation, the step size At is restricted by stability requirements. However, in this present model, At is only restricted by the degree and amount of nonlinearities present in the system. This restriction is far less stringent. C. Practical Limitations and the Cures It was necessary to solve several difficulties in the present simulation study in order to provide practical solutions. We shall describe them item by item. 1) The electron energy distribution is the quantity of major interest. The difference between the magnitudes of its maximum and minimum non-zero values is of the 9 . order 10 . This is due to high density in the low-energy region (thermal energy) and the relatively low density 113 in the high- energy region (near the source energy) . It is unrealistic to hope to produce reasonable numbers of particles in the low-density (high-energy) region 9 and at the same time have 10 times as many in the low- energy region. This problem is resolved by combining the following techniques: a) Russian Roulette and splitting. b) Weight factors which produce fractional particles. c) Logarithmic energy intervals. The number of particles in each energy class interval AE . is directly propor- tional to f(E.)*AE. (f(E.): number density at E.). By using logarithmic energy intervals, we have small intervals at low energy region and large intervals at high energy so that the product f(E.)*AE. assumes a practical magnitude. 2) The total number of simulating particles (machine particles) N T is limited by the storage requirement and the economy of the calculation (computing time is propor- tional to N ) . We have to set an upper limit (Nmax) on the total number of simulating particles in any time cycle.. During each time cycle of duration At, N source particles are introduced into the system, and N secon- dary particles are produced (by ionization) . At the same time N particles escape from the system, and N , particles are lost due to recombination. There are two conditions to be satisfied at all times: 114 a) N + N + N ,FID<60l,NREZtlll,EZI< UI.NTPM Or02 DIMENSION FNEI55 I , FNREZIU I, NRZ 1 1 111 ,CNE< 35 I , CNREZI 11 1 ,TNE (551 000 3 COMMON 0, E(5100),Z15100I, UI5100) , I AS< 500 I , IRRI 8501 ,NRR , IES( 600) X,WT(5100I,SWF 0004 COMMON /AREA1 / NIZN , NOZN, NAS, 01 ZNI5 I ,N SE, NMET, NEXC, ETAS , XWNME'T", WNEXC, WNIZN," ELS, UNAS ,WNSE 0C05 COMMON /AREAZ/ CI ,MUO, E I ,EI S, SF, ANEAM,PQP0,PI ,E0,NI0N, J E , JZ tNASWl , 1NSEW ,AKT, AIKT , AMI , AMC1 , EDI , PDMUC, P0?6,K1 , K2~~t K3. MN2.EDZ,ED3 0C06 COM MON /AREA3/ XSEL , XELC , XRCB, XMET. XEXC ,XIQN.XT0L,KELC 0007 " COMMON /AREA4/ CONST, 0EL2, ELT,IE, Oil, DIZ, 013 ,DELI 0008 COMMON /A5/ EG <_55 I , E0EI33I, EBDI55I 0009 COMMON >A6/ ESEI50I , NUP ,NDNR,NDN, WSETSOI 0010 REAL MUO. NION.ND C STATEMENT FUNCTION FOR EVALUATING MAXWELL IAN AND RELATEO FCN? 0011 F M1(£,FNEM(-EXPI-E»AIKT I »SQRT Tall — RNJfUZlWU&SSllI S0»P0P0*SR*1.0E*16 SF-l.O AITG* CMUPDP /t AKT»PI> OE1-0.6T6 DE2-0.355 DE3-0.5 DII-1./0.016 D12-1./DE2 D13-1./0EJ ■ _ AMC1- SQRTf2.0*2.7l83*AIKtl*E02 DO 11 1-6, 36 EBL «( I-5l*DE2 -4.0 EMV- EBL - 0.1685 EBD,( I I- 2.0**EBL EGU-1]^ 2.0**EMV 11 EDEl : ll-EBDUI*I2.C**DE2 -1.01 JE1-JE+1 DO 12 1-37. JE1 EBL" 7.J>05_*Jl-36I*DE3 _ EMV- EBL - 0.1405 EBD( 1J^ 2.0**EBL_ EG(I-1»- 2.0**EMV 12 EOEI I>- EBD1 1 >»(2.0»»0E3 - 1.01 DO "10 I "It 5 EDE( I )» 0.016 _ EG(I )» EDEl 11*1 -0.008 10 EBD< II-_EDE(II*( I-l> EDEI36I-29.69 E GI48I- 173 7.482 _ IF (NTAPE.NE.U GO TO 2 READ(10» EAVG. DAVG, TME,_NTW, NT, HEADI10I Z, WT, ANI0S.SF3 NS«NSW*K2 _ WGHT-1.0 REWIND 10 NTP, NION. DOEG. NSW. SWF, E,U NT«NT*1 2 CONTINUE IF (NTAPE.NE.2I READI9 I EAVG, RE AD (9 I Z, WT, .. NS_ -HSW«K2 WGHT-1.0 REWIND. 9 NT-NT*1 5 CONIJNUE_ _ IF (NTAPE.NE.Ol GO TO 5 DAVG, THE, NTW, NT, NTP, ANI0S.SF3 NION, DDEG, NSW, SWF, E, U GO TO 122 128 FORTRAN IV. G LEVEL 10 MAIN DATE • 7*126 05 tiUSi 0139 0140 0141 0142 0143 0144 0145 0146 0147 0148 0149 0150 0151 0152 0153 0154 0155 0156 0157 0158 0159 0160 0161 0162 0163 0164 0165 0166 0167 0168 0169 0170 0171 0172 0173 0174 0175 C176 0177 0178 0179 0180 0181 0182 0183 0184 0185 0186 NION- SQRTIS0*10/AITGI AN£AM-N10N»0.33333/OEN CALL S0RCE2INTP, EAVG, NTH, PM, FID» A L 'A I QA( EAYG) OAVG-AL DOEG-DAVG«ND/SORT(EAVG) SF2-1.0 IF (NION.NE.O) 1EHMJM/H1W. SN-S0*AL»ND*SF2/!ETV*SQRTIEAVGI I NSH- A NSW »0.5 IF (SN.GT.ANSWt NS.NSM»K2 NSW-SN +0.5 SWF-NSW/SN SF-SMF SF1-SF2 SF3-SF2 CALL SOURCEINS.NTP) TflE-QiO. ANIOS-NION HGHT-1.0 , 122 CONTINUE C GENER ATE SOURCE _P AR T J CL_ES_ J 0_ INITIALIZE THE PROGRAM ANEAM-NION*0. 333 33/OEN WRITE16.70) 70 FORMAT!' 'IS "PROGRAM INITIALIZATION AND THE INITIAL PARAMETERS* ,//> H R1TE16.71) A 1TG. NI ON, NTP . AL.6AVG 71 FORMAT*' ', 'AITG -'.E14.6, 3X,'NI0N - • , E 14. 6, 3X , *NTP -•,112,31, l 'OCENlEOi -'. F14.6.3K. 'EAV G- '.F14.4.n . 72 FORMAT!' ' , 8X , • OCOL • , 1 3X , • OCEN' , 13 X , • ENERGY' , I IX , • I • , 14X , ' U • , 10X, . „,;. X'MEIGHT'./I 75 FORMAT!* *,F14.6,3X, F14.6, 3X, F14.5.3X, F12. 4.3X, F12. 4,3X, El 3. 41 M R1TEJ6.72I . RIP-RAN3Z<0> C WITHIN THE OUTER I LOOP THE~ I TH PARTICLE HISTORY r r.FOMFT Rir. rhutinf ENTERING here 7777 IF IIEiII.LE.O.OI.OR.IEin.GT.EO*1.5M GO TO 9008 ... TCEN-DDEG/.EJX __ 7778 IF UNIZN.GT.OI.ANO.iI.GT.NTPM TCEN- DIZN!N0ZN*1I Ul I ) -U!II»IUSN ... ...__ ._. LFG2-0 ETA.G1-1.Q IF !E( D.GT.EDlt ETAGi — 1.0 SQEje _iQBJ!ElIJ> ... _ . _ DVA-EF*TCEN*EMASSI*0. 5 VZA"Sfl£*E.TV*Ul !)_..♦ DV.A . VYS-E! I»*ETV*! 1.0-U! I»*U! 1 1 »*ETV nr.FN- S0RT!V?A»«2+VYSI»TCEN NFLG-0 qnn TF IFdl .LF. O.OI £0 TO _ 9Q08 IF (NFLG .EQ.OI GO TO S09 DVA-EF*TCEN*EMASSI*0. 5 129 FORTRAN iv G LEVEL It , MAIN DATE • 72124 OTmTO 0187 VZA-50E*ETV*Um ♦ DVA ' 0188 809 DZ- VZA»TCEN 0189 ENEW" E(II ♦ EMDZ 0190 IF ( ENSW.LE.O.OI GO TO 9008 0191 SQEP-SQRT(ENEM) 0192 UN EW» SQE«U > 0199 "" DCflL- ALIUALD 0200 IF (NFL G.EQ.lt GO TO 901 0201 ~ IFTmod ( I , I M I . EO. » lrf"RTTTf677TT DCOL.DCEN.EU ) ,Z( I ) ,U( ! ITVTfTI 0202 DC 0L2-ABS(AL0G(1 < -RI l)*ALD 0203 " 0T0T-0T0T ♦DC0L»0C0L2 02 04 ET0T-ET0T*E1II 0205 IF (N0PN.LT.2I 615 TT5 5oT 02 06 TCEN2»TCEN 0207 DCEN2-0CEN 0208 EIO^EIU 0209 UIO» U(I»*IUSN 0210 ZIQ-Ztll 0211 ENMO-ENEW 0212 UNW0-UNEM_ 0213 OZO-DZ " 0214 901 NFLG"! 0215 902 IF (OCEN. LT.DCOU GO TO 905 0216 TCEN-TCEN - 0COL»TC£N/0CEN 0217 DZC-DZ*DCOL/DCEN 0218 E(IL"E(_I»^EF»_DZC 0219 zm-zii > *ozc 0220 DC EN- OCEN -DCOL 0221 NC0L«NC0L*1 0222 I F (MO0INCO LiMC».EQ. 0» RIP-RAN3ZI0) 0223 IF (MOO(NCOLfMC).NE.OI RIP-l.O-RIP 0224 CALL COLPROU, TCEN, NTP, RIP, NUST, &9008, C9045I 0225 IF((E( I).LE.6.0).ANO.(LFG2.EQ.l) t GO TO 9008 0226 GO TO 900 _ 0227 905 E(I»»ENEW 0228 U- UNEW _ 0229 Z.OR.(Z(n.GT._0_n GO_TO 9008 0231 j-INT(Z(I)*DZIVi*l 0232 IF (J.GT.JZI J»JZ 0233 K» INOEXE(E(III 02 34 __ IF (K.GT.JEI K- JE 0235 IF (K3.E0. II GO TO 904 0236 ETAG2 -1.0 0237 IF lEIII.GT.EDir ETAG2— 1.6 0238 IF <(ETAGl*ETAG2I.GT.0.0J GO _T0 _?04 0239 IF (ETAG2.GT.0.0I GO TO 903 130 FORTRAN IV G LEVEL 18 MAIN DATE - 7212* 05/31/62 0240 0241 NUP-NUP*1 ESEINUPI-EUI 0242 0243 WTU I-WTMI/K3 WSEINUPI-WTUI 0244 0245 GO TO 904 903 IFULFG2.EQ.il. OR. IN0PN.LT.2M NONR-NDNR*! 0246 0247 IFKLFG2.EQ.il. OR.INDPN.LT. 211 NDN-NDN+1 IRR(NDNRt-I 0248 0249 IF (M0D(NDN t K3).NE.0l GO TO UTIII«MTm*K3 9045 0250 0251 904 NEZtK.JI- NEZ(K.J) ♦ MTdt EKZ(K,JI-EK{(KiJI*Em*WTm 0252 0253 9045 IFK LFG2.EQ.il. OR. (NDPN.LT.2M LFG2-1 GO TO 9009 0254 0255 EK1-EI0 um«uio 0256 0257 ZK l-ZIO TCEN-TCEN2 0258 02 59 DCEN-DCEN2 0C0L-0C0L2 0260 0261 ENEW-ENWO UNEW*UNWO 02 62 0263 DZ-DZO IF KUSN.LT.O) GO TO 900 0264 02 65 GO TO 902 9008 NES-NES^l 0266 02 67 IF (E( II.LT.ED1I HTK I-WTK1/K3 IF IWTlII.GT.il WTm-1.0 0268 0269 HNES>WNES ♦WT(I) IESINESI-I 0270 0271 IF (EIII.GT.O.O) ETEP-ETEP*E(I»»MTIII 9009 I« 1*1 02 72 0273 LFG2-0 IF (I.LE.NTPI GO TO 7777 0274 0275 IF K I.GT.NTP).AND.(NIZN.GT.O) t IF INOZN.LT.OI NIN-1 NOZN* NOZN-1 02 76 0277 02 78 0279 IF KNIZN.GT.Ot . AND. (NI N. EQ.O II . _ JLLZW»H_NJjt_N__»Q.5 . NWS-HNES +0.5 NASW-WNAS *0.5 GO TO 7778 02 80 0281 WNAS-MNAS/KELC WNSE-WNSE/KELC 0282 0283 0284 0285 NSEH-WNSE +0.5 _QQ L9JL_ N-LrJE. - _ .._.. 00 191 N-l.JZ EE(HI-EE 0297 22 IF (NK3.GE.JZt IZ-JZ 0298 DO 21 L-ltIZ 0299 NREZ(LI-NREZ«L» ♦ M3*WSE(Ni 0300 21 EZ(Lt«EZ(L> ♦ ESE t N>«N3*MSE(N» 0301 NK3-NK3 -JZ*H3 0302 IF (N K3.LT.JZ> IZ-NK3 0303 N3-1 0304 IF (NK3.GT.0) GO TO 22 0305 23 CONTINUE 0306 2* CONTINUE 0307 DO 192 M-l, JZ 0308 DO 1 915 N-l.JE 0309 NREZ(MI-NREZIM) ♦ NEZIN.MI 0310 1915 EZ(H>-EZ(HI ^EKZIN,*) 0311 192 CONTINUE 0312 NTPN-NTP*NIZN-NES- NAS 0313 NTPW»0 0314 E TOL-0. 0315 00 14 M-l.JE 0316 ETOL» ETOL »EE(H> 0317 14 NTPW-NTPW *NE(M» 0318 EWTOL-ET OL 0319 NTPT-NTPN 0320 DAVGO- OAVG»NO _ 0321 DAVS-(DAVGO*1000.»**NTVP 322 05AVS»0SAVS* DAVS 0323 DTAVS-OTAVS *OAVS 324 _ DI5G"1.0/05AVS _ 0325 ANTPO - 1.0/NTPW 0326 N5TP- N5TP ♦ NTPW »0.5 0327 ~NTP5«NTP5* NTPN*DAVS 0328 N TPO-NTP 0329 EWAVG»ENTOL*ANTPO 0330 NSO»NS 0331 NSW0>NSO*WN2 *0. 5 0332 NTP- NTP+NIZN 0333 NTPP*NTP 0334 DAVG'0T0T»0.5/NTP 0335 EAVG- ETOT/NTP 0336 DE G-DO EG/ETV _ ____ 0337 TME"TME*DEG 0338 T HS-THS*PEG C IF SOURCE PARTICLES SOINSK NES fHEN ELIMINATE AS KANV RESIOUE PARTICLES AYTH 0339 NTPS»NTPW* NWS +NASW 0340 NTU^NTPM ♦. 5 0341 PCE N1 - AL0A1 EAVG I _ 0342 DDEG*DAVG*ND/SQRT( EAVGI 0343 BF-1 .0 _ 0344 IF HWNES+MNASI.GT.O.OI 8F« ( HN I ZN+SN I / (HNES+WNA S» 132 FORTRAN IV G LEVEL IB »UIN OATE • 72126 05/31/02 0345 SF 2*1.0 0346 u IF (AN10S.NE.0I Sf 2-NTPS/ANIOS , 0347 ESftVG-0.0 0348 IFIHNES.NE.O.OI ESAVG- ETEP/MNES 0349 EASAVG-0.0 0350 IF IWNAS.NE.Ot EASAVG-ETAS/WNAS | 0351 SHFO-SWF 0352 CN5S-CN5S* SN»DAVS»WGHT 0353 CSF1- CSF1 ♦< ETEP*ETAS*40. *MNIZN l/TNS 03 54 SFi»CSFl«DI 5 G/tS0»E0> 0355 CSF2- CSF2 ♦ SF2*0AVS 0356 IF ((H0D(NT,5).EQ.0».AN0. (NT.NE.Ot) SF3-CSF2»Q 136 0357 WGHT-SF1/SF3 0358 SN«S0»0DEG*SF3/ETV 0359 IF (NTP.GT. 49201 ANSW-2 0360 NSWANSW »0.5 0361 IF (SN.GT.ANSMI NSM-SN +0.5 03 62 SWF-NSH/SN 0363 NS-NSM*K2 0364 SF-SHF 0365 NMEM>MNMET *0.5 0366 NSXW-HNEXC +0.5 0367 M5VLE-M5VLE* HVALUE*0AVS 0368 CETEP-CETEP «■ ETEP»DAVS 0369 CETAS-CETAS *ETAS*0AVS 0370 C N 5IZ-C N 5 U»WN It N»PAVS - 0371 CNES*CNES*MNES*DAVS 0372 CNAS.CNAS»WNAS»DAVS 0373 CNMET«CNNET*WNNET*OAVS 0374 CNEXC-CNEXt>MNEXC»DAVS 0375 CESAVG-CESAVG* ESAVG*0AVS 0376 _ CEASCCEASG*EASAVG»DAVS 0377 CLOSS- CL0SS*IWNIZN-WNES»*0AVS/ CNREZIM! ♦ NREZ(MI*DAVS 16 CEZ(M|- CEZ(M) ♦ EZ(MI*OAVS 00 18 M-l.JZ FNREZ(M)«NREZ(M)*ANTPO 18 EZ(N)« EZ(M(/ NREZ(M) SUMNE-0.0 DO 193 M-l.JE SUMNE- SUMNE ♦ NE(M» TNEIMI- TNEIMI ♦ NEIMI/EOEIMI 193 FNE(MI-NE«M)*ANTPO/EDE EAV3-0. DO 1755 M-l, JE EAV3- EG(MI»FNE(MI*EDEIMI *EAV3 1755 CNEIMI- TNE(M»«AN5TP CNION-CNUPDP*2.0*AIKT/SQRT(AKT*PI» ANION-0.0 ANI02-0.0 DO 20 H-l.JE ANI02- ANI02 * FMHEGINI ,CNE HF«t ANIZVG*ANSVG»/( ANES+ANAS > IF (MODCNT.NPMI.NE.OI N ION-ANIOS GO TO 1790 IF (IWF.GT. 0.0051. AND. IWF.LT. 200.11 NION-ANIOS*WF**NBF 1790 CONTINUE DAVT«DAVG*ND 13U FORTRAN IV G LEVEL 18 MAIN DATE - 72126 09/31/02 0448 0449 0450 0451 0452 0453 0454 0455 0456 0457 0458 0459 0460 0461 0462 0463 0464 0465 0466 0467 0468 0469 0470 0471 0472 0473 0474 0475 04 76 0477. 0478 0479 0480 0481 0482 0483 0484 ALPHA-UN IZN/(NTPM*DAVT I ANEAM-NIQN»0.33333/DEN 1785 FORMAT CO' ,' ANION- '.E 13.4, 3X,'DNI ON-', El 3.4. 3X.-TNI0N-' ,E 13.4, 3X, X'ALPHA-'.F12.5.4X.»LEAK-'.E13.4./I : URITEI6. 17851 ANION. 0NI02, TNI02. ALPHA* ALOSS WRITEI6.171I EO. PDPO. D. OETZ. PEG. NDPN 171 FORMAT!'!' ,'GASt HE •, 2X, • SOURCE ENERGYi ', F5. 0, 3X ,' PRESSURE i • ,F4. 0, X 'T0RR'.3X.'SLAB TH1C K NESS> j .E*. 1 . ' CM',' DXt',F4.1,« CH' ,4 X,'PEL T X=',E12.5,3X,'»PROC-',I2,/ I WRITEI6.18 0I NT. THE. NS.NSU 183 FORMAT!' •, "NUMBER OF TIME CYCLES-' . 16, 3X , 'TIME -', E12.5,' SEC, 13X. 'GENERATED SOURCE PART -', 15. 3X, 'MGTED * OF SOURCE PART.-', 115 I L8___l _B AVT __ 0CEN1.EAVG.EWAVG. EAV3 1805 FORMAT!' ','DAVG -'.F12.5.3X, 'OCEN -• ,F 12. 5, 3X, • OVERALL EAVG-'. 1F10.5.3X.' WGTEO E AVG-' .F10. 5.3X. 'EAVG C HECK-' .F10. 5 ./I WRITE(6,181I 181 FORMAT!' '. IX. 'NUMBER 01 STR1 BUTION' ,2X. 'FLUX DENS 1TV .4X. ' AVERAGE IE ENERGY', 3X, 'CUMULATED # D1STR' , 3X, 'ENERGY INTERVAL' ,2X, 'CUM FLUX XDEN' ,4X.'«PEN-PEVIATI0N './I . CESUM-0.0 DO 281 M-l. JE . CESUM-CESUM ♦ CNEIM)*EDE!M) IF INE1MI.GT.0I PEAY6«eG<'H-EE> / < EGI w yk i _. C C I U _» 111 ? B1 URITE16.182I H.FNE I H I . FLX1 . EGI Ml .CNE! HI . EDE ! H I .FLX2 . PNAVG 182 FORMAT!' • , 13 ,E12. 5.5X.E12.4, 5X, F12.6.8X, E12.5.6X, F12.6.4X.E12 X.5.4X.E12.4I DET-ODEG/ETV MRITE16.183I NTPP. N5TPN.N TPW .PET. DAVGO. NCOL 183 FORMAT!' ',// ,3X, • TOTAL • OF PARTICLES PROCESSEO -• , I5.5X, 'CUMULAT _ 1FP i OF PARTICLES -' . I10.3X . ' WGTED TOTL • OF PAR T ICLES-' .F12. 3 .//. X' NEW DT-',E12.5,4X,'0LD DAVG-' ,F12. 5.4X, • • OF COLLISIONS-' , 15 , //I yRlTE16.187l 187 FORMAT!' • ,1X, 'ZONE NUMBER'. 5X, 'AVERAGE ENERGY' ,5X, 'NUMBER IF FLFCTR0NS'.3X. 'CUMULATIVE » PEN ( I I ' . 3X . 'CUMU. A VG ENERGY'. 3X. l'ZONE » OENSITY' ,/l .DO 282 M-l. -J-Z ANZ-CNREZ!MI«AN5TP*PI5G PC7-CF7IHI/CNREZ!HI 282 WRITE(6,184I M, EZIMI, NREZ(M), ANZ. PCZ,FNREZ(M» 135 FORTRAN IV G LEVEL 18 mmr DATE - 72124 05/31/0r 04t85 0486 0487 0488 0489 0490 0491 0492 0493 0494 0495 0496 04 97 0498 0499 0500 0501 0502 0503 0504 0505 05 06 0507 0508 184 FBKHrm '. **♦ I S, SX. FU.6,M,FU.t, fT, F U.8. M.FU.l. T I, — 1 F12.5) - " WRITE(6,183I NION,WN!ZN,NIZN, ANIZV6,WNES. NES.ANES, MN A S.N AS, XANAS t NTP, Kl t K2, HF.BF 185 FORMATC'O", 'NUMBER OF POSITIVE ION-' ,E 14.6,// , • NUMBER OF ION I Z AT 1I0NS 0CCURE0-',F9. 3 , 5X, 'ACTUAL !■', 15, 3X, 'CUMULATIVE AVC»» ,F10. 3, x//,' » of ESCAPED PARTICLES-', F9. 3,8y,. ' AC T UALI- ' . I S, *X, ' CUMULA TI VE X AVG«',F9.2,//.' # OF ABSORPEP PART1C-' ,E12. 4 , 3X, 'ACTUAL t-',I 5,3X X, 'CUMULATIVE AVG- ,F9.2 ,//, • TOTAL I OF PARTICLES IN THE NEW CVC1T XE -',15,//,' K1-',I2, 3X,'K2-',12, 4X,» BALANCE FACTOR WF-", XF10.5,3X,'BF2«', F10.5,/» WRITE(6,188> E01. NPM, NO, KELC Ta e FORMATS ENERGY PlvlPlNG LINE FOR RR PLAV E - ' . F 10.4,4X, ' EV E RYMZ X,'TH_ CYCLE CHA NGE N »' ,//,' TIME INTERVAL MU LTIPLY FACTOR-' ,F9. 4, X3X, 'ELASTIC XSEC / BY', I 2,/ I NOR - N DN-NDN/K3 WRITEI6.186I NSO, NSW, ANSVG,NOR,WNSE,WNMET, ANMET.WNEXC, ANEXCNUP, XWGHT 186 FORMAT!' • • • SOURCE PARTICLES GENERATEO OLP -',110,//, 1' SOURCE PARTICLES (WEIGHTED! IN THE NEW CYCLE -',U0,3X. 'CUMULAT XIVE AVG-',F10.3,// 1 ,' * OF PARTICLES KILLED OFF IN RR PROC-' , IIP,//, 1' « OF SUPER ELASTIC COLLI SION-' ,F10.4, // , 1 ' • OF EXCITATIONS TO META STABLE LEVEL-' X.F10.4, 3X, 'CUMULATIVE AVG-' , F9.3, I //,' i OF EXCITATIO NS TO OTHER LEVELS-' ,F 1Q. 4, X4X, 'CUMULATIVE AVG-' , F9.2 ,// , ' • OF PARTICLES CROSS OVER TO Hl'SH XE NERGY REGION NUP- ' , I 5.4X, ' OVERALL WT-' ,E12.4, / I WRITE(6,195I EF, EOP, POPO ,SF1, SF2, SF3 195 FORMAT I' E-F1ELD-', F9.2,'V/CM', 4X, ' E/P-' , F9. 2, 4X, 'P/PO -',F7.1, X3X, 'SF1-', E12.4,3X, • SF2-' , E12. 4, 3X, 'SF3-', E12.4./I MRITE(6,170I S0,WVC,SWF0,WVALUE,WVLE1, WVLE2,ESAVG, AESAVG,EASAVG, I-', XAEASG 170 FORMAT! ' SOURCE RATE -',E11.3,« PER SEC ' ,4X, 'W-VALUEI LOCAL XF10.3, 3X, 'SOURCE REP. FACT-' ,613.4,//, • W-VALUE-', XF10.3.4X,'CUM AVG WV1-' ,F9. 3,3X, 'CUM AVG WV2-',F9.3, 1 //,' AVERAGE ENERGY OF ESCAPEO PARTICLES-' , F12. 5, 3X, 'CUMULAT XIVE AVG-'.F11.4,//,' AVERAGE ENERGY OF ABSORPEP P A RT 1 CLES-' ,F12. 5, X3X, 'CUMULATIVE AVG-' ,F11.4,/ I CKFNE-SUMNE+ANTPO WRITE(6,189>CKFNE, CESUM, TNION, TNI02 189 F ORMATC ' , ' THE SN » CHECK ■' , F12. 4, 3X, ' THE CN « CHECK-' ,F12.4,//, 1« CALCULATED N* USE F-' ,E14. 6,3X, 'CALC. N* USE CN -',E14.6,3X, »AR X TIF CN+-'.E14.6,//) IF ((M00(NT,5».NE.0».0R. ( MODI NT, 10 ». EQ.O M GO TO 3 WRITSI9) EAVG, OAVG, TME , NT W, NT, NTP, NION.OPEG, NSW, SWF, E, JJ WRITEI9I Z, WT ,ANIOS, SF3 EWFUE 9 REWIND 9 3 CONTINUE IF (MOO(NT,10).NE.O) GO TO 4 WRITEdOl E AVG, DAVG. TM E, NTW, NT, NTP , NION, PPEG, NSW, SWF, E, U WRITEUOI Z, WT ,ANIOS, SF3 136 FORTRAN IV G LEVEL 18 MAIN DATE • 72126 OS/31/02 0309 0510 0511 0512 0513 051* 0515 0516 0517 0518 0519 0520. 0521 0522 0523 0524 0525 0526- ... 0527 0526 0529 0530 0531 Q5J2 — 0533 0534 0535 0536 0537 0538 0539 0540 0541 0542 0543 054it. 0545 0546 0547 0548. 0549 0.5.50 0551 0552 0553 0554 0555 0556 0557 0558 C559 0560 ENOFILE 10 RFMlNO 10 * CONTINUE 1F n| - wg - Q> 6666 CONTINUE NT- MT ♦ 1 ftp Tfl 6666 00 56 MM 1 500 nWW tMi-0.0 IES(N»"0 JABJJUafl 56 IASIMI-0 p p si M-1.SO HSE(NI"0.0 11 F^F' MI "°- Q I«l norHT l Al WF ™* OUTPUT QUANTITIE S DO 52 M»li JE FFIHI»0. 52 NE(MI- 0.0 nn S3 BslaJl NREZ(M>« 0.0 C7IMI-0.0 DO 53 N«lfJE FKPM.MfO.O 53 NEUNtM»»0.0 137 FORTRAN IV 6 LEVEL 18 MAIN DATE - 72126 05/31/02 0561 NIZN-0 0562 NOZN-0 0563 NIN-0 0560 NES-0 0565 NSE-0 0566 NCOL-0 0567 NAS«0 0568 NASW1-0 0569 OTOT-0.0 0570 NNET-0 0571 NEXC-0 0572 WNES-0.0 0573 WNAS-0.0 0574 WNEXC-0.0 0575 WNIZN-0.0 0576 WNMET-0.0 0577 WNSE-0.0 0578 ETAS-0.0 0579 ETEP-6.6 0580 ETOT-0.0 0581 NW-0 0582 NUP-0 05 83 NON-0 05 84 NONR-0 0585 INIV-INIV* 1100 0586 CALL RN3INZIINIVI 0587 WRITEI6.72I 0588 IF (NT.LE.NTCI GO TO 7777 0589 9999 STOP 0590 DEBUG SUBCHK 0591 END 138 FORTRAN IV 6 CEVET Tf" TOTETF DAT! - TlI2i 09/31/02 0001 0002 0003 0004 0005 0006 0007 oooe 0009 0010 0011 function indexeiei COM MON /AREA*/ CONST. 0EL2. ELT,tE, Oil. 012. 013 ,DEL1 IF (E.LT. 0*081 GO TO IF (E.GT. 158. 1335 1 GO 11 TO 13 12 tNOEXE- 5*(AL0G(E»*C0NST ♦4.01*012 RETURN 11 INDEXE-l ♦ E*DI1 RETURN 13 INOEXE-36 ♦ CALOG6 r*>n\n oou? ODOS <*«■> ■>'• : 7 00 J« 0010 r -11 -'■12 t-'^n CM* 0015 CU 0017 c .: 1 H 0019 0020 r 21 0022 0023 0025 CC26 0027 OC26 0C29 0030 Of 31 0032 0033 00 34 0035 0036 0037 0038 0039 0040 CC41 0042 0043 0044 0045 0046 0C47 0048 SUGKOUTlNfc Sr)PCE2(NS0t EAVG, NTWtPM, FIOI "LAL Nl(55),Fin(50> ,WI(55» Common I), r(5ie>l,2(5lCCI, UtSlOOl . 1 AS1 5CC » , IftUI 850 I t N»H , I* S< 6 * : I X,«T( 51001 ,SWF COMMON /AKEA2/ C1,MU0,EI,FIS,SF,ANEAM,P0I>C,PI ,frf , NION , JF , JZ ,NA*MJ , INSEM ,AKT, AIKT ,AMl , AMC 1 ,E01 , PONUO, Pf)26,K1, «?. Kit N^iEDNEOl COMMON /A5/ EG(55I, E0H55I, E80(55> FMt(t)« tXP(-E*A!KT)*SOt-T(E»«AHl FM2(M»AA/f r : M3(F»-EX(M-»:*AIKn*SURT(E» »AMC1 KEAL MIJO. NION <>04 FURMATd «,F12.6,3X, 13, It, E12.4, 2X, F 12. 6, 2X, H2. 5 , 3X, H 2 4 I 605 F(I«MAT(« •,• Jl-« ,I3,3X,«E01»»,F9.3 t ?X , »NSO • , 1 *, 3X ♦ «K1»» , 1M.3X, 'K2*«, I«,//,' INITIAL DI STR 1 WITHIN* ,//, • #0F PTS • , 7X, • I NOCX l',5X, •« OfcNSITY* ,4X,»AVG ENERGY* ,6X , • OFL-t • ,*X, • *F IGHT* ,/ I 6"6 FORMAT! • •,«NSl« , ,I5,*X,«NS2" , ,I5t3X,»NSO- , .l5,3X, , NTW»« . 15,3X, X«NTW*EXP(-PM| Jl«lNDEXE(Eni-0.5l J2»Jl*l ISUM-0 HRIT5I6tf05l 00 1"? 1-l.Jl ONI -F10II »*EOE(I )-NTw/Kl Mil l-DNI wl( I 1-1.0 IF IFIDm. (ONI.GT. (ONI, (DM, f UNI, Jl, EDI, ►'SO, Kl, K2 112 IF IF IF IF r,T, LT, ,LT, 112 100 202 20: LE.C.O) GC TO 220. I Nil 11*200 200.1 WI(l l-DNI/ 200. l?l NKII-12 12) MI(II«0NI/12.0 NRITE(o,604) NKII, I, FIDdl.EGdl, EDEIII.Mllll I SUM- I SUM* NI (I I JSUM-0 DO 20^ I-J2, JF. DNI «FIOm*EDE(I »*NTW*K2 NI (I MOM Mid 1-1.0 IF (FlOd I.Lc.O.Ot GO TO 202 IF (ONI.GT. 203. I NKU-200 IF (ONI.GT. 200. I Ml ( ! l-ONI/200. IF (0NI.LT.12) NIdl-12 IF (ONI. LT. 121 Wl (I l-ONI/12.0 WRITE(6,604I NKII, I, FlOd ),EG(II, EDFdl, Midi JSUMx jsum ♦ NKII NSUM«ISUM*JSUM NT-2-ISUM-K1 +JSUM-WN2 *0.5 MRITEI6.606I ISUM, JSUM, NSdM, NTW, NTM2 NTW-NTW2 IE-i EAG-0.0 00 300 I -I, JE lUO FORTRAN IV li LEVEL Id SDRCE2 OATF - 72126 05/31/02 00*o IF (Nim.LT.ll GO TO 300 C0 5C JI-NKII ♦!£ -1 0051 00 301 J-IE, JI, 2 0052 R!afUN3Z(* FBOU 1 ♦ FOE( I»*R2 0056 wT( JI'HM I 1 0057 -T(JM I'UII [) 0C5P F.AG -EAG ♦ F(J) +EIJUI 005*WTU) 00 72 too CONTINUE 0073 EWAVG-ETW/NTW2 007* HRITE(6,607» EAVG, EMAVG 0075 QDZ«D*0. 25 0076 I»l 0077 101 RZ«RAN3Z(0I 0078 RU-RAN3ZI0) 0079 zm»o*Rz 0080 UU)-2.0*RU -1.0 0081 I- I* 1 0082 102 Z(II«D»RI 0083 um»i.o- 2.0*RU 008* l«_l..»l ... . . . 0085 103 HI »«0*(1.0-R2» 0086 UII)« 2.0*RZ-1." 0087 I- I *l 0088 10* ZII)»D*I1.0-RZI 0089 um» l.-2.*RZ Q090 1* It 1 ... . . 0091 IF(( 12 -I-3».GE.0» GO TO 101 0092 IF (I2.LT.I) GO TO 105 0093 RZ»RAN3Z(0) 0094 RU-RAN3ZI0) 0095 IF(( 12 -M.EQ.2I GO TO 102 0096 . . . ... . I£iJ-12_ill*ifl.li_ GO TO Lfl3.__ 0097 IF(< 12 -II.EQ.OI GO TO 104 0098 105 RJrTUBJ* 0099 END EAG»FAG-EU»n 1U1 FORTRAN IV G LEVEL 18 SOURCE DATE » 72126 05/31/02 0001 SUBROUTINE SOUPCE I NSt NTP I 0002 COMMON 0, E<5100). ZI51001, UI5100I . I A S( 500 > , IRR < 950 » ,NRR , I E S( 600 J X,WT<51CP|,SWF 00m COMMON /AREA2/ CI . MUO , E I , E I S , SF , ANEAM.PDPO ,01 , FO , NION , JE , JZ ,NA SWl , 1NSEW ,AKT, AIKT .AMI ,AMC1,ED1,P0MU0, PD26.K1, K2 . K3, WN2,feD2,E03 0004 KEAL MUO, NION C SOURCE PARTICLE PARAMETERS GENERATING PROCEDURE C II: THE INITIAL INDEX VALUE FOR THE SOURCE PARTICLE 0005 I*NTP*1 0006 DO 100 K-li NSt 2 0007 E(II"EO 0008 RZ«RAN3Z(0I 0009 RU-RAN3ZICI 0C10 Z(M«D*RZ 0011 U(II«RU*RU-1.C 0012 HT( I l»WN2/SWF 0013 IF UNS-IO.LE.OI GO TO 100 0014 E(I«1)-EC 0015 Z I I ♦ 1 » » RZ»D 0016 U< 1*1 1— Ul I I 0C17 taTI I*ll*rtN2/SMF 0018 I»I*2 0019 100 CONTINUE 0020 NTP«NTP*NS_ 0021 105 RETURN 0022 END 1U2 FORTRAN IV G LEVEL 18 MAIN DATE • 72126 05/M/02 0001 0002 0003 0004 0005 0006 0007 0008 0009 0010 0011 0012 0013 0014 0015 0016 0017 0018 0019 0020 0021 0022 0023 00 24 C WHENEVER THERE ARE" MORE HOLES THAN" SOOTCT P ARTICLE'S, FTLL IN THE EXTRA HOLES W C WITH PARTICLES AT THE END SUBROUTINE FILHOLtNTP, NES, NASI COMMON 0, EI5100),Z(5100I, UI5100I , I AS< 500 I , I RRI 850 I ,NRR , I ESI 600 I X,wTI5l00l,SWF 100 105 300 200 K-NTP IF (NES.Lt.OI DO 100 1*1 • J-IESI II E(JI- EIKI Z(J>* ZIK) U( Jl« U( Kl HTIJI'WTIKI K*K-1 NTP-NTP-NES IF INAS.LE.OI K*NTP DO 300 I»l, J-IASI II e< j i -e ( k i Z( J>*Z(K> UIJI'UIKI WT(JI>MT(KI K*K-l_ NTP-NTP-NAS" RETURN END GO NES tO 105 RETURN NAS 1U3 FORTRAN IV G LEVEL 18 RRPROC DATE » 72126 05/31/02 0001 SUBROUTINE RRPROC INTP.K3I '"" 0C02 COMMON D, E ( 5 100 I , I ( 5 100 ( , U(5100» , US < 500 I , IRP ( 850 > ,NRR , I ES ( 600 » X,WT(5100I,SWF 0003 COMMON /A6/ ESEI50I , NUP .NDNR.NON, WSEI50I 0004 NR-NDNR 0005 K-NTP 0io 6 IF (K3.LE.ll RETURN 0007 100 CONTINUE 0008 I«IRRINR» 0C09 E 0010 UIII'UIKl 0011 Z(I)«Z(K) 0012 WT(It>Wt(K) 0013 K«K-l 0014 200 NR-NR-1 0015 IF (NR.GT.OI GO TO 100 0016 NTP«K*l 0017 300 RETURN 0013 END Ikk FORTRAN IV G LEVEL 18 ALDA DATE - 72126 05/31/02 00C1 00 2 00C3 0004 0005 0006 0007 0008 0009 0010 0C11 0012 CC13 0014 0015 0016 0017 0018 0019 0020 0021 0022 0023 0024 0025 0026 0027 0028 0029 0030 0031 0032 0033 0034 0035 00 36 0037 FUNCTION ALOA(E) COMMON /AREA2/ CI , MUO, E I , F I S , SF , ANEAM, PDPO, P I , EO , NI ON , JF , J* ,NA $W 1 , 1NSEW iAKT, AIKT .AMI , AMC1 , EDI , PDMUO, P026.KI, K2 , K3. WN2.ED2.ED3 COMMON /ARFA3/ X SEL , XELC ,XPC B, XMET , XEXC , X ION, XTOL ,KELC REAL MUO, NION WELCM.O/KELC IF (E.LE.5.T XELC-23.6*PDP0 IF< (F.GT.5. I. AND. (E.LE.19.BI I XEL C-23. *PDPO« SQfi T ( 5 . /E I IF (E.GT. 19.81 XELC-EXP(-0.02*E»*16.76*P0P0 GO TO 5001 - 19.801**2 - 1 9 ;i oT7 rei ♦ DIMENSION PX(4t FM3(EI-EXP(-E*AIKT»*SQRT(E)»AMC1 REAL MUOi NION C DECIDE WHICH KINO OF COLLISION XTI -ALDAIEdll if (em.LE.eo3 » go to 5006 C ENERGY REGION 1 19.0EV TO I KEV 4000 PMET* XMET*XTI PEXC* XEXC*XTI _ PION* XION*XTI PELC«XFLC*XTI PX<1 1* PEXC PX<2)» PEXC *PMET PX(3»* PXI2I +PION IP»1 DO 45 IX«1,3 45 IF(RIP.GT.PX( IX) I 1P-IX+1 GO TO (4015.4013*4017,40111, 4011 Rl-RAN3Z(0t U(I>« 2.0*R1 -1.0 RETURN _ C INELASTIC COLLISION PROCESS EXCITE C EXCITATION TO METASTABLE LEVEL 4013 WNMET* WNMET ♦MTU I Ed)- Ell l - 19*80 IF (E(II.LE.O.O) RETURN1 NMET»NMET»1 IP TO META STABLE STATE RI 1> RAN3Z(0I Udl - Rl l ♦K11-— - 1.0 RETURN C EXCITION TO O THER LEVELS 4015 WNEXC -WNEXC +WT ( I I Em» mt -2it»o IF <€< II. LE. 0.01 NEX C»NEXC*1 _ RU* RAN3Z(0I U(II - Rll ♦*!! ■ RETURN IPJUlIIIOfl P«QtSSS RETURNl l*Q 0040 0041 0042 .M43_ 0044 4017 WNIZN- WNIZN* WTdt _EJ_I1*_ Ed) - £J IF (EUI.LE.O.OI RETURNl N1ZN«NIZN»1 PRODUCE SEC ONORV ELECTIONS NOZN-NOZK ♦! 1U6 FORTRAN IV G LEVEL 18 COLPRO DATE ■ 72124 05/31/02 00*5 00*6 0C*7 00*8 00*9 0050 0051 005? 0053 005* C055 0056 0C57 0058 0059 0060 0061 0062 0063 006* 0065 0066 0067 0068 0069 0070 0071 00 72 0073 007* 0075 0076 00.77 0078 0079 0080 0081 0062 008* 008* 0085 0086 0087 0088 0089 0090 0091 0092 0093 C IONIZATION PRODUCED PART" "lCL"BS~ A"lWAVS~PuT fS't'HF ENO~ OfHlfSIDUES DIZN(N01N»- OCEN . _ . RI-RAN3ZC0I UIII»KI*R|-l.fl ..._ . NNTP «NTP *NIZN E(NNTF.i- £111*0,1 Z»Z< I I U(NNTP)— U(I» E( I) ' E( I »*0.9 WT(NNTP)- MT(I) IF ((E(NNTP) .GT.EDD.OR. ( (E( I I»1.11UE1 I.LT.F01 ) I GO TO 4011 NDN«N0N*l NDNR"NDNP«-1 WT(NNTP>"WTU»*K3 .. IRRINDNR)-NTP*NIZN *018 CONTINUE RETURN C ENERGY _RE_G_JJJJL I LQlRl«FMIEl7"~ FR2«FH3(£2I IF (RE.LE.FR2I RS-RE . 1F__CP2.IE. FP 1 ) g feHI IF < DELE-Cl»(Em-ENNI»(1.0-C0SSTt*ELS»KELC Ell I- EUI -DELE um-ui 50125 IF( ELS. LE. 0.01 BJLtUJLN. C RECOMBINATION PROCESS U ( I I - 2.0*RAN3Z(0» -1.0 1U7 FORTRAN IV G LEVEL IS 009* 5013 NAS-NAS+1 0095 IASINASI-l 0C96 WNAS- WNAS*WT(l) 0097 IF (E(I).GT.O.OI 0098 PETURN2 C SUPERELASTIC COLLISIC 0099 5015 E ( I ) » EI 1 1 ♦ E I S 0100 WNSE-MNSE ♦WTI II 0101 NSE«NSE+1 01C2 ftP« RAN3ZI0) C1C3 um- RP *RP -1.0 0104 RETURN 0105 ENO COLPRO DATE ■ T2126 03/31/02 ETAS«ETAS*E(II*WT(I I PROCESS 1U8 FOR TRAN lv a li :vul in c MASTE" P L TICL r S I c DEPIN1TI r 0: l)IA«E c. i 1 : INI T I c JZ : Ti'Tai c IAS: AfcP bn: IJ 1ST c V: JUNE c ASSAYS 6 r NTR: T IT c NT; TOTA c c" mi MIMH NIZN: Nil occi kE 1*101 0002 DIME or 3 COMrt CO 04 X,wJ( ~ comm" XwNME. 000? COMM INSEw 0006 COMM 0CG7 COMM OOOfi COMM 0009 COMM 00 10 REAL OX 11 c STATEMEN FMU 0012 DATA MAIN DATE - 72125 C1/2B/29 0013 .?CJ.i. 0015 0016 0017 DATA DATA 1.0/7 OATA KilliOAM OF MONTE CARLO SIMULATION OF NONLINEAR ALPHA «>LASM»S (ALOHA NfMJCED PLASMA FL tCfRON" (tNTRBY D ISt«" I ftUT ION STUDY ONS AND TFKMINOLOGY TER OF THE CYLINOEft. THE LFNGTH OF THE FINITF. REGION 4L ENERGY OF THE SOURCE PARTICLES SC(NSM SOURCF RATE NUMBER OF ZONES JEJ TOTAL NUMM.R OF ENERGY RANGES AY FOK STORING THE INDEX OF ABSORPEO PARTICLES NAS« T HE NO. OF THFM ANCF TO ROUNDPY ' OCOL~: 0"lST ANCf'TO COLLfSloKT 0"CE"N~: OtsrANCE TO C£N INTERVAL . Z. li STOPING ENERGY. DISTANCE. AND OIRECTION OF PARTICLES At NUMHER OF PARTICLES TO BE PROCESSED EACH CYCLF L NUMBER OF TIME CYCLES f-K OF SOURCE PARTICLES PERJJNIT T IM?*DI STANCE OR SOURCF PATE MBER OF IONIZATION PROCFlS OCCURRED tACH TIME CYCLE AL NF(55(, NEZt 55, 10 » , F I i 1 1 » , DX ( 3 ) , EE< 55 I , CEZl 111 tEKZf 19 ,FlC(6d» ,NRFZ(U).EZK11»,NTPW NSIPN FNf(55», FNREZ(ll), NRZ I (1 1 ) , CNE( 55 » t CNRE/ I 1 1 ) , TNF « 55 » ON 0. F<510C >,Z(5100I, U(5100» , I AS ( 500 ) , IRR( «550 I .NRP , IF S< 600 » 51001, SWF ON /ARFA1/ NIZN, NOZN, NAS, 01 ZN( 50C j ,NSE% NMET, NEXC, ETAS , T. wNfXC, WNIZN, ELS,MNAS ,WNSE ON /AREA2/ Cl,MU0,fcI,EIS,SF,AN£AM,POPC,PI,EO,NION,JE,JZ,NASKl, ,AKT, AIKT ,AM1 .AMC1.E01, PDMUU, PD26.KI, K2 , K3, WN2.E02.ED3 ON /AREAV XSEL,XELC,XRCB,XMET,XEXC,XION,XTOL,KELC ON /AftFA4/^ CONST, DEL2, ELT.IE, Oi l, DI2, DI3 .DELI ON /A5/ EG{55»,~ EDEI55I, EB0<55» ON /A6/ ESEI5CI , NUP ,NDNR,NON, WSE(50I MUO, NION.ND T FUNCTION FOR EVALU ATING MAXWELL IAN jAND RELATED FCNS E,FNtM»*EXP<-b»AIKTf*SQRTfE)*FNEM TM S, NT, DTOT, E TO T, NI N, N I ZNT,NAST/0. 0,1,0.0. 0.0 ,0,0,0/ w5VLE,CN5IZ,CN5S,CNES,CNAS,CNMET,CNEXC,CESAVG,CEASG/9*0."*7~ oio SUS CNF/ 55»0.C/, CEZ/U+0.3/, CNR E Z/l 1*0. 0/.TNE/55»0./,Elt?/5S0»Q D5AVS/0.b/ NEZ/ 550*0. /,NE/55 *0./. EZ/11»0./.NRSZ/ 11*0,/ , EE/S5*0.0/ READ! 5,1 I NTAPE P EA D(5,15 II Kl, K2, NTC, NTP,NTW,NPM, KELC ,JZ, 1M,NCUM, NDPN,NTMP 0014 0055 0056 QQ5 7 Cl»2.'*0. 010549/4. 003873 ETV* SORT (2* 0*1. 60203/9. 108 1 *l.0E*08 VTf« 1.0/«FTV*ETVI CALL UNDERZI 'OFF* > P026=PnPU*26.0 MU? = < 1.0F-C9I/ETV r.rtlJ0 = 1. 0F-09 CMJPDP=CMU0*PDP1 IF ( POPr-.CT. 10. C > P9MU0*C"ijPDP/ETV CLOSS=". r NUP*0 _"i!iN3j .... _ CMUP0P«>CHIJQ*1Q,Q 0058 CC59 0060 OQfil 0062 0063 C064 -O&Ai.. 0066 -0C67 0068 006.9 00 70 _OQJJ_ 0072 0073 0074 OPTS 00 76 0077 0078 0079 0080 nofli 0082 0083 0084 0085 NDNP ■? INO»0 NAS = NES*0 NSb-J... NG7N*0 NIN»0 LK = iMv^^o^pgg NMET=0 NEXC=0 ETAS = r>.C ETEP-0.0 NTP 5-0 NW»0 N5TP*0 DTAVS=0." _C_SJ=2i0^0 CNI02=0.0 CETAS=0.0 f.FTFPsP.P NCOL-0 . HNME T'Q .O. WNSE«0.0 WNEXC^O.0 WNES»0.0 MNAS^fKP WNIZN'O. JE«4d _. IE*JE ,DEN»3.-5AE*L6»RDPQ_ DETZ=D/JZ CONST=1.0/ALOG(2.0I _ED2»0. 447513 . E03«19.8 PM«Fn?»A!KT 150 FUPTBAN IV G LEVEL lfl MAIN OATE - 72125 01/2H/29 C086 oca'1 00R8 QC89 0090 0C91 ecu C093 0095 SO 96 Of 9 7 009A TC99 0100 0101 cn? 0103 010* 0105 0106 C1C7 0108 01C<5 0110 0111 0112 0113 0114 0115 0116 r 1 1 7 0118 C119 C12T 0121 0122 0123 0124 0125 0126 0127 0128 AMl-7.r*AIKWSQRT(PI»AKT) CALL f>M3INf(?"«»<»T632rn — SO-POP^* SK*1.0E*16 AITG- TMUPOP /< AKT»PI I DEl-O.Olb Pt2« - .?55 DE3*0.3 11 0129 0130 0131 13? 0133 9P * 0135 0136 011*1. /(J. 016 DI2-1./OE2 cn«i./ne3 AMC1- SOkT(2.r'.*2.71fl3*AIKT>»ED2 00 11 1=6, 36 EBL = < I-i>)*D£2 -4.0 EmWTBL" - r .l685 FBO( I (- ?.0*»EBL rG(I-l)« 2.0**EMV FDF( 1 )rfn0( I )*(2.0**0E2 -I. 01 JE1*JE*1 00 12 1=37, JE1 _ _ ""pr'l'b" 7.005 ♦(I-36I»DE3 FMV= Cf'L - C.140 5 FHDl I )= 2.0»*EBL FT,( 1-1 >» 2.0-**E*V 12 EDE= EBP< I »*<2.0**DE3 00 10 1=1, 5 FOE( I 1= 0.016 EG< I )= EDL( I )+I -C.008 10 EBD( I t= EDt( I »*( I — 1 » E0E( 3bl=29.6V EG(4ri»» 1737.482 _IF_ (NTAPJE.NE.il GO PtADIlM EAVG, DAVG* PEADI10I Z, wT NS*NSW*K2 _J*EWINO 10 NT«NT*1 2 CONTINUE 1.01 TO TME, NTW, NT, NTP, NION, DDEG, NSW, SWF, E,U ANI0S.SF3 IF (NTAPE.NE.2l GO TO 5 REA019 ) EAVG, DAVG, TME, NTW, NT, NTP, NION, DOEG, NSW, SWF, E, U i)EAD(9 | Z, WT, ANI0StSF3 NS=NSW*K2 _ _ REWIND 9 NT*NT*1 CONTINUE IF _iNTA_PE.NE.O) GO _T0__ 122 NION- SQBT OAVG=AL OOEG = DAVG»ND/SOP.T I EAVG) 0137 0138 SF2=l.O IF (NION.NE.Ot SF2-NTW/ N10N 151 FORTRAN IV G LEVEL 18 .. MAIN DATE "72125 01/28/29 0139 SN»S0»AL*ND»SF2/(ETV«SQRTt£AVGI ) 0140 NSW«ANSW*0.5 0141 .... IF ISN.GT.ANSw.J . N5W«SN_.*0j5_ 0142 NS»NSW*K2 0143 SKF?.NS-W/SN 0144 SF-SWF 0145 ANIQS'NKN 0146 CALL SOURCE(NS,NTP> 0147 SF1-SF2 014H SF3-SF2 0149 TME=O.C 0150 122 CONTINUE _D1A1 SF=S*F GENERATE SOUkCE PARTICLES TO INITIALIZE THE PROGRAM 0152 ... . ANFAHx,mion*0.333?3/DE l N .,_._■_ 0153 WRITE(6,70) 0154 70 FORMATCl't'PROGRAM INI T 1 AL 12 AT IO N AND THE INJTIAJ. PARAMETERS',//) 015^ WMITE<6,71I AITG, NtON, NTP, ALiEAVG ill 56 71 FORMAT!' '. 'AITG = '.E14.6, 3X.'NIQN - ' , E 1*.6, 3X, 'NTP -MlZ .n, l'DCbN(EO) =', F14.6.3X, 'EAVG» ',F14.4,/> 0157 72 FQKMATC • , 8X , 'OCOL' , 13*, •D£EJ4A»ilX J .' ENERGY* , 11X, • I • , 14X t «U' , 10X, X • WEIGHT', / I C15H 75 FORMATC »,F14,6,3X, F14.6, Mi F14.5,3_X, F12. 4 ,3X ,Fl-2,4, 3X , Fl 3. 4) 0159 WRITE(6,72I .0160. K1P = »AN3Z(0J 0161 1=1 I WITHIN THE OUTER 1 LOOP THE I. T_H. PART ICLf HISTORY C GEOMETRIC ROUTINE ENTERING HERE 0162 _ .....77J.7...1F.i.LLlll.LE,0.0UQA*iE(IlAiiIxlfl.U GQ LO 9008 0163 DCEN*OOEG*SORT(E(I I) ri64 7776 IF MNI/N.GT.OI.AND.U.GT.NTPH DCEN» PI ZN< NOZNM I 0165 NFLG=0 0166 0JJJ=U11J*1USN _.. . 0167 LFG2-0 ETAG1-UQ_ . .. 0169 ni 70 moo IF !E< II.GT.ED1I ETAGl—l.O IF 1(711 l.l F.O.OI.nR.(ZtII.GT.DI.QR.(E(I I.LE.O.OI 1 GO TO 9008 0171 01 7? J»INT(Z(I)*DZIV» *1 RI.RAN17C0) 0173 r»t 74 ALD«ALDA(E( It I 0COL*ALD*ABS(ALQGIRI 1 1 0175 0176 IF (NFLG.EO.l) GO TO 809 / IFIMnnil.lNI.EQ.OI WRITEI6.75I t>COL. DCEN.EI I 1 .Zf I 1 .Uf I 1 .HTI I 1 0177 DC0L2*ABS( ALOGI 1. -Rll l«>AL0 01 7ft DTOT-DTOT t0C0L»DCQL2 0179 ET0T-ET0T*E(I I ni«n IF 1N0P N .LT.2I GO. 10 609 0181 DCEN2-DCEN nift? Fin-El II 0183 UI0« UIII*IUSN 0184 ZIDx/lll 0185 809 NFLG-i 0186 900 IF IDCF N .LT^DCQL) GO 10 902 152 FORTRAN IV G LEVEL 18 MAIN DATE ■ 72125 01/2H/24 _01B7 901 OCFN'OCFN- DCOL 0188 ZIII-ZUI ♦ OCOL*UUI ' '. 0189 NC0L«NC0L»1 0190 IF <«.OD(NCOLiMC).EO.0l RIP"RAN3Z(0I 0191 IF (MOD(NCOLtMCI.NE.OI R IP* 1 ■ O-MP 0192 CALL COLPROUt DCENt NTP, R'T>t NUST, C9608, H9C09I 1 9? Gn TO ?oe o , 0194 902 2(11" Z < I » ♦OCFN*inTT 0195 IF < INOfcXfc(E< I II 0199 IF (K.GT.JE) K-JE 0200 "' IF"(K3.E0.n GO "TO "SO* 02"! ETAG2-1.0 _ 020? IF (E( I I.GT.E01I ETAG2--1.0 0203 , IF < (ETAG1*ETAG2).GT.0.C| GO TO 9C4 C2G4 IF (ETAG2.GT.0.0) GO TO 903 0205 NUP»NUP*1_ 0206 " estiN"n>i«E(n " ' •--•—-- 0207 WT(I l*WTU I/K3 020A HSF(NUP)«wT( I I C2C9 GC TO 90* 0210 903 IF< < LFG2.EQ.il. OR. JNDPN.LT. 211 NDNR-NONR^l C 2 1 1 _J FKLF G? . b . 1 I . R .(NDPN .LT.2II IRR(NONRI«I 0212 "" NDN«N0N*1 0213 IF ( MnD(NON,K3).NE.OI GO TO 9045 0214 WY(I l=WT(I l»K3 0215 904 NEZ(K,JI» NEZIK.JI ♦ WT(M 0216 Ek76 11 IF (NK3.GS.JZI IZ-JZ 0257 00 21 L-l.IZ 023% NRFZ(L1«NREZ(LI ♦ M3»wSF(N> _ r2"j9 ?1 ^Z ( L » -E7 I L » ♦ FSE(N»*«*3*WSE(NI 0260 NK3-NK* -JZ*M3 0261 IF (NK3.LT.JZl IZ-NK3 0262 M3«l 0263 IF (NK3.GT.0) GO TO 22 C26f. ,?3 CONT INUE 0265 24 CONTINUE f'66 DO V->2 M»i, JT 0267 DU 1915 N»1,JE 0?6d N»F.Zm»NHFZ(MI ♦ NEZ(N.M) 0269 1915 EZ =CZ ( M I ♦FK/JN.MI D2.7.Q LK CU-NTI^UX . ._ C271 NTPN*\*NIZN-NFS-NAS 0272 NTP.J-0.0 027* FTUL-0.0 0274 L)Q 14 V=1,JE 0275 ETUL» FTCL ♦C-FJMI . 02_7a. 14 NTPd.NTPW »NF(H1 , 0277 EWT(JL*ETCL 0278 ..N.TPT»NTPN 0279 UAVr,n= 04V0»N0 0280 DAVS'IOAVGOMOOO. )**NTMP 0281 D5AVS=JSAVS* OAVS A2.a2 D_L5_G = U1ZQ1AY5 0283 ANTPO = 1.0/NTPW 0234 N5TP- N5TP ♦ NTPW *C.5 0285 DTAVS-DTAVS tOAVS 0286 NTP5*NTP5» NTPN*OAVS 0287 NTPO=NTP Q2B8 cwAvr.«F^ Tni*ANTPn . 0289 NSO-NS C29C NSWQ*NS0*wN2 *0.5 0291 NTP«NTP*NIZN G292 NTPP*NTP 15U FORTRAN IV G LEVEL 18 HMN DATE • 72125 01/28/29 0293 "0294" 0295 0296 0297 0298 02 99 C3 0C 03->l 302 C3C3 0304 "0 3 05 C3 06 0307 306 *■> "3 0310 0311 0312 0313 0314 0M5 '■316 -"317 02 18 on 9 T320 0321 0??? 0323 0324 ^3?5 0326 0327 3 2 8 0329 0330 0331 0332 33J 0334 0335 0336 0337 0338 -0J3JL 0340 0341 342 0343 DAVG"OTOT*0.5/NTP EAvfc- eToT>nT^ DEG-ODFG/ETV TMF"TME*DEG TMS«TMS*OEG C IF SOURCE PARTICLES SOINSX NES THEN ELIMINAtC AS MANY RESIDUE PARTICLFS ' NTPS«NTP* ♦ NWS »NA5W DCFN1- ATOM' EAVGI NTW-NTP* *,5 OOEG«DAVG*ND/SQRT< EAVGI . SF2-1.0 IF (ANIOS.NE.O) SF2-NTPS/ANI0S ESAVG-0.0 IF(WNES.KE.O.O) EfAvTT-~l"«P/WNP5 FASAVG-ETAS/WNAS BF»(WNIZN*SNI/IWNES*WNASI EASAVG-O.C IF ( WNAS.ImE.O) irtFD»SwF PF«1.0 IF ( (wNtS + WNASKGT.OjiOl CN5S*CN5S+ SN*DAVS CSF2- CSF2 +(ETEP*ETAS*40.*WNIZN»/TMS IF(IM00(NT,5t.E0.C). ANO.INT.NE.CN SFI-CSF2 /(SO*EC*S.I CSF1- CSf-1 ♦ SF2*DAVS IF I jK? > 0361 NAST« NAST ♦NWS C INTRODUCE THE SOURCE PARTI CLES FOR THE NEXT CYCLE C EVFBY FIFTH CYCLE RENORMALIZE NION ANO OUTPUT QUANTITIES ARE PRINTED C RfcNUkMALI/ATICN PROCESS OF NION i?6? AN5TP* 1.0/N5TP 0363 ANTP5-05AVS/NTP5 ji&A N5TPN*N5TP r>3oi... J2a._i6-._M-l i JZ 0166 CNHEZ(M)* CN«EZ(*> ♦ NREZ(M»*DAVS 0367 16 C E _ ( ^ I = 0C7(M) ♦ EZ(MI*OAVS 0J6J 00 l« M-l.JZ 0?69 FNREZ , FNEl M) >*EDE (M » 0386 20 CONTINUE 0387 . DNION*ANICN*CNIJN 0388 DNI02=ANin2*CNION _3fl9 ___Q_AL__S . ._.._. 0390 IF ( ALCS-.LT.C.OI CJD-C.0 0391 _■_ IF 1 A_S(C_DI.GT.11Q.*SQ) CJD-ii__^SQ 03<»2 TNI0N*ARS(S0*CJ0)/DNI0N 0393 ... _ _INID2-ABS(S0*CJDI/DNIC2 156 FORTRAN IV G LEVEL IB MAIN OATI ■ 71125 Cl/21/2* 03 94 TNI02-SQRT(TNIQ2I . ~0~393 '' TNION-SORTJTNlONI ' '. C396 CNI02-CN102 ♦ TNI02*OAVS 0397 ANr0S»CNI02*0I3G 0^98 wF-1.0 0399 IF (ANES4ANAS.NE.0.0I WF-I AN I Z VG*ANSV£ 1 1 \ ANIS*AN AS I 04 00 IF (MOO!NT .NPM}.NE.O) G O TO 1790 C4C1 NION-ANIOS 040? IF ! (wF.GT. 0.005 LAND. (WF.LT. 200.1) N10N-ANI0S»WF**N*F 0*03 1790 CONTINUE 0404 OAVT«DAVG»ND 0405 ALPHA«wNlZN/lNTPW*DAVT> 04C6 ANEA M »NICN»Q. 33 333 /PEN ^_^ 0407 1795 FORMAT! '0* , 'ANION-' ,TT3.4,3X, 'ONION- * ,1 IT. 4 t TxT 7 t"NI0N« , ",ETi^*", i*t X'ALPHA-' , F 12. 5, 4X,' LEAK-', El 3. 4,/ I lv)8 ' hP lTfc(6,1785 ) ANIONi ON 1 02, TNI02, ALPHA, ALOSS 0409 wklTt (6,171 > EO, POPO, 0, DETZ, OEG, NOPN 0410 171 FORMATCl' ,'GAS: ME ', ?X ,' SOURCE ENERGY J ', F5.0, 3X ,' PRESSURE! ', F4.*, X_ •jn°« , t3X 1 _« SLAB THICKNESS: '.F^Jj,' CM',' 0XI',F4.1 L » CW ,4X, 'OfcLT "~ X-"' ,E12.5,3X,'»PROC«' , 12,/ I C411 WRITE!6,18GI NT, TME, NS.NSW 0412 140 FORMAT!' ', 'NUMBER OF TIME CYCLES- ', t*>, 3X, • TIME •', *12.5,' SEC, 1»X, 'GENERATED SOURCE PART -', 15, 3X, 'HGTtO « OF SOURCE PART.-', 115 I 0413 wRITE<6tl905l OAVT, DCEN1 , EAV G, Ew AVG, EAV3 -41^ 1H:5 FORMATC ','OAVG - • , F12.5, 3X , •DC'EN - • ,~F"1T. 5", ?X, • OVERALL PAVG- • , 1F10.5.3X,' WGTEO E AVG-' ,F10. 5,3X, 'EAVG CMFCK-' ,F 10. 5 ,/ I r'415 WRITE!*, ISl) 0416 141 FORMAT!' ', IX, 'NUMBER D I STR IBUT ION • ,2X . 'FLUX DENS ITV • ,4X, • AVFRAGF It ENERGY' ,?X, 'CUMULATED « DISTrt ', 3X, • ENERGY INTER VAL ' ,?X, 'CUM FLUX XDEN' ,4X, '«DEN-DEVIATION%/l ' 04 17 - - - - (• ESUMs0#0 OMk 00 281 M-l, JE n^i<; CESUM-CGSUM ♦ CNE ! Ml *EOE (Ml 0420 TF (Nc(M).GT.OI OEAVG-EGI MI-EE( Ml /NE . ANO. ! QE AVG.GT .Q. 0) I .OR. ( (M. EO. JE I . ANO. ( OEAVG.LT.0.0 I I X) DNAVG--0EAVG*FNE(M|*0.5/EDE(M| 0429 IF ( ! .E12.4,/ | 0t5C WRITEJ6,195) EF, EDP, PDPO ,SF1, SF2, SF3 0451 195 FORMAT C E-FIELD"', F9.2,'V/CM'» 4X, • t/P«' , F9. 2, 4X, • P/PO .',F7,1, X?X, »SF1 = ', !!12. 4, -*X, 'SF2«',F12.4,3X, 'SF3-', El?. 4, /I .C!»,>,», »AR XTIF CN*-',E14.6,//I •' IF nHOD.lQ.OH 60 TO 1 WRITE!*! EAVGi BAVOt Tll» NTM, NT, NTfi "NTONifiOIC, NSW.fwf, f, u 0*57 0*5 8 0*59 0460 0461 0462 0463 0464 0466 0466 0467 0468 0469 C470 0471 0472 047? 04 74 0475 047t 0477 04 78 0479 C4 8C 0481 04H2 0-4 83 0484 0485 04 86 0487 04«8 0489 049C 04Q1 04 92 0493 0494 0495 04 96 0497 0499" 0499 0500 0501 0502 05P3 WRITE<9> I, MT ,AN10S, ENOFILE 9 REWIND 9 3 CONTINUE SF3 IF (M00_(NT,10I.NE.0I GO TO _*_ MITeTfCI EAVG, DAVG, TM£' ( NTw, NT, NTR, NION, OOEG, NtW.SNF, E, U WRITE(IO) I, WT .ANIOS, $FS_ ENOFILE" 10 PEWlNO IC 0504 0505 0506 0507 4 CONTINUE IF (MOOtNT,NCUP).NE.0> N5TP-0 NTP5-0 D5AVS«C.C NAST-0 NIZNTaC K5VLE*C.r CLOSS«0.0 CSFl-i)." CNI02«0.0 __CN5IZ-0.0 CN5SiC*. > ' CNM?T»0.0 CNfcxc*c .c CESAVG«0,0 CEASG-0.0 00 54 M«1,JZ CNREZ(MI«0.0 54 CEZ(M)-C.7 OU 55 M«1,JE TNE(MJ=0. 55 CNFMI*r." IF MUD(NT i 10).NE._0> '"'UfAVS«0."0' CfcTEP*C.O rETAS-0.0 CNtrSO.n CNAS-0.0 6666 _CONTJ NUE NT« NT" Vl 00 56 **1,500 OIZNJMJxC.r JRRIM)*? IES(M|«0 56 I A S I ^ > = " _ N«l,50 GO TO 6666 GO TO 6666 00 51 __WSf " 0.0 C511 --00.. _51...M-lt.il 0512 NRE2(M|« 0.0 csn __iH)-0_C- ___ 0514 00 53 N-l.JE 0515 n_iii_.«j»_j_a 051b 53 NE2IN_(INIV) Ob** IF (NT.Lh.NTO GO TC 7777 C543 9°99 STflP 0544 DEBUG SliUCHK C545 END i6o VITA Benjamin Shaw-hu Wang was born in Inner Mongolia, China, on March 25, 1937. He is the son of Mi. C. H. Wang, the former Chinese Army General now the Congressman of the Republic of China. He received the B.S. degree in Electrical Engineering from the National Taiwan University of Taipei, Taiwan, China in 1960, he then served as an ensign in Chinese Naval Academy for one and a half years. He came to the United States in September, 1962 and enrolled in the University of Missouri at Rolla, Missouri. After he received the M.S. degree in Electrical Engineering in 1964, be joined the Westi nghouse Research Laboratories and served as a research engineer for three and a half years. In September 1967, he came to the University of Illinois at Urbana-Champaign to undertake studies for a Ph.D. in Computer Science Department. He has been a research assis- tant and teaching assistant in the Department of Computer Science and a research assistant in the Nuclear Engineering program. He is the member of Kappa Mu Epsilon and Pi Mu Epsilon the two, national honor societies of mathematics. LIOGRAPHIC DATA ET 1. Report No. UIUCDCS-R-72-519 3. Recipient's Accession No. 5. Report Date May, 1972 Itle and Subtitle Monte Carlo Simulation of Nonlinear Radiation Induced Plasmas uthor(s) Benjamin Shaw-hu Wang 8. Performing Organization Rept. fllUCDCS-R-72-519 irforming Organization Name and Address Department of Computer Science University of Illinois at Urb ana-Champaign Urbana, Illinois 6l801 10. Project/Task/Work Unit No. 11. Contract /Grant No. Sponsoring Organization Name and Address Department of Computer Science University of Illinois at Urb ana-Champaign Urb ana, Illinois 618OI 13. Type of Report & Period Covered Ph.D. Thesis 14. iupplementary Notes ibstracts A Monte Carlo simulation model for radiation induced plasmas with nonlinear properties has been developed, which employs a piecewise linearized predict-correct iterative scheme to resolve nonlinear characteristics. Several variance reduction techniques have been developed and incorporated into the model which include the antithetic variates technique. The model has been applied to the determination of electron energy distribution functions in a noble gas plasma created by a-particle irradiation. The model can be applied to other types of complex plasma problems with nonlinear properties. ey Words and Document Analysis. 17a. Descriptors /•lonte Carlo Simulations, Monte Carlo Calculation of electron-energy distribution functions, Monte Carlo method for radiation induced plasmas, electron energy distribution functions in noble gas plasmas, toite Carlo Techniques for plasma problems, •lonte Carlo Simulation model for nonlinear plasmas. identifiers /Open-Ended Terms =• OSATI Field/Group /ailability Statement "elease unlimited "►■ms-35 ( 10-70) 19. Security Class (This Report) UNCLASSIFIED. 20. Security Class (This Page UNCLASSIFIED 21. No, of Pages 22. Price USCOMM-DC 40329-P7 1 JUL 2 6 WW • < I V w I I H1H ■ I ■ i ! ■ HI