UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Particle simulations of a one dimensional plasma Alfheim, Arne 1983

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.

Item Metadata

Download

Media
831-UBC_1983_A6_7 A43.pdf [ 4.54MB ]
Metadata
JSON: 831-1.0085213.json
JSON-LD: 831-1.0085213-ld.json
RDF/XML (Pretty): 831-1.0085213-rdf.xml
RDF/JSON: 831-1.0085213-rdf.json
Turtle: 831-1.0085213-turtle.txt
N-Triples: 831-1.0085213-rdf-ntriples.txt
Original Record: 831-1.0085213-source.json
Full Text
831-1.0085213-fulltext.txt
Citation
831-1.0085213.ris

Full Text

PARTICLE SIMULATIONS OF A ONE DIMENSIONAL PLASMA by Arne A l f h e i m B . S c , C o n c o r d i a U n i v e r s i t y , 1980 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE i n THE FACULTY OF GRADUATE STUDIES Department of P h y s i c s Plasma P h y s i c s Group We a c c e p t t h i s t h e s i s as c o n f o r m i n g t o the r e q u i r e d s t a n d a r d s THE UNIVERSITY OF BRITISH COLUMBIA OCTOBER, 1983 © Arne A l f h e i m , 1983 In presenting t h i s thesis i n p a r t i a l f u l f i l m e n t of the requirements for an advanced degree at the University of B r i t i s h Columbia, I agree that the Library s h a l l make i t f r e e l y a v a i l a b l e for reference and study. I further agree that permission for extensive copying of t h i s t h e s i s f o r scholarly purposes may be granted by the head of my department or by h i s or her representatives. I t i s understood that copying or p u b l i c a t i o n of t h i s t h e s i s for f i n a n c i a l gain s h a l l not be allowed without my written permission. Department of The University of B r i t i s h Columbia 1956 Main Mall Vancouver, Canada V6T 1Y3 Date DE-6 (3/81) i i ABSTRACT A P a r t i c l e In C e l l (PIC) code has been d e v e l o p e d t h a t s i m u l a t e s a one d i m e n s i o n a l plasma. The plasma model was e l e c t r o s t a t i c w i t h f i x e d i o n s . The code was used t o s i m u l a t e Landau damping and t h e two stream i n s t a b i l i t y . R e s u l t s o b t a i n e d showed good agreement w i t h the t h e o r e t i c a l b e h a v i o u r f o r b o t h phenomenon. C o n s e r v a t i o n of energy was a l s o checked. The t o t a l energy was found t o d e v i a t e o n l y s l i g h t l y so t h a t energy c o n s e r v a t i o n was a c h i e v e d . i i i TABLE OF CONTENTS 1. INTRODUCTION 1 2. THEORY 4 I n t r o d u c t i o n 4 Debye S c r e e n i n g 4 Plasma O s c i l l a t i o n s 5 Landau Damping 7 Two Stream I n s t a b i l i t y 11 3. DESCRIPTION OF PIC CODES 15 I n t r o d u c t i o n 15 The Model 19 N o t a t i o n 23 The Scheme 24 Some Important S c a l i n g Parameters 31 4. RESULTS 34 LANDAU DAMPING 34 Energy C o n s i d e r a t i o n s 48 R e v e r s i b i l i t y 53 TWO STREAM INSTABILITY 57 5. CONCLUSIONS 68 REFERENCES 71 APPENDIX A . .72 APPENDIX B . . . 74 LIST OF TABLES i v T a b l e I . T h e o r e t i c a l Damping Rates And F r e q u e n c i e s For Landau Damping. 11 T a b l e I I . T h e o r e t i c a l Growth Rates And F r e q u e n c i e s For The Two Stream I n s t a b i l i t y 14 T a b l e I I I . Landau Damping Rates O b t a i n e d From S i m u l a t i o n s 44 T a b l e IV. Sample E n e r g i e s And D e v i a t i o n s 52 T a b l e V. Computed Growth Rates F or The Two Stream I n s t a b i l i t y 64 V LIST OF FIGURES F i g u r e 1 Scheme For P a r t i c l e Mover 25 F i g u r e 2 Time H i s t o r y Of The P e r t u r b a t i o n i n K Space D e n s i t y P e r t u r b e d 39 F i g u r e 3 Time H i s t o r y Of The P e r t u r b a t i o n i n K Space Thermal V e l o c i t y P e r t u r b e d 40 F i g u r e 4 L e a s t Squares F i t To Peaks I n F i g u r e 2 41 F i g u r e 5 y/u T h e o r e t i c a l And Computed 46 F i g u r e 6 u/u T h e o r e t i c a l And Computed 47 F i g u r e 7 Energy H i s t o r i e s ( E l e c t r o s t a t i c and K i n e t i c ) Thermal V e l o c i t y P e r t u r b e d 49 F i g u r e 8 Energy H i s t o r i e s ( E l e c t r o s t a t i c and K i n e t i c ) D e n s i t y P e r t u r b e d 50 F i g u r e 9 T o t a l Energy 51 F i g u r e 10 T o t a l Energy For D i f f e r e n t S i z e d Time Steps....54 F i g u r e 11 Time H i s t o r y Of The P e r t u r b a t i o n i n K Space V e l o c i t i e s Reversed 56 F i g u r e 12 Phase Space P l o t For Two Stream I n s t a b i l i t y kv =0.97 59 F i g u r e 13 Phase Space P l o t For Two Stream I n s t a b i l i t y kv =1.1 61 F i g u r e 14 K i n e t i c and E l e c t r o s t a t i c Energy For The Two Stream I n s t a b i l i t y 62 F i g u r e 15 T o t a l Energy For The Two Stream I n s t a b i l i t y . . . . 63 F i g u r e 16 K i n e t i c and E l e c t r o s t a t i c Energy For The The Two Stream V e l o c i t y . V e l o c i t i e s Reversed...66 F i g u r e 17 T o t a l Energy For The Two Stream I n s t a b i l i t y . . . . 67 v i ACKNOWLEDGEMENTS I would l i k e t o thank my s u p e r v i s o r Dr. B a r n a r d f o r h i s h e l p and guidance t h r o u g h o u t t h i s work. For answering my many q u e s t i o n s c o n c e r n i n g A s s e m b l e r , I would a l s o l i k e t o thank L u i z da S i l v a . I g r a t e f u l l y acknowledge t h e f i n a n c i a l s u p p o r t g i v e n by the N a t i o n a l S c i e n c e and E n g i n e e r i n g R e s e a r c h C o u n c i l and the U n i v e r s i t y of B r i t i s h C o lumbia. 1 CHAPTER 1 INTRODUCTION To i n v e s t i g a t e the p h y s i c s of plasmas, t h e r e a r e t h r e e avenues of approach. The two t r a d i t i o n a l methods a r e t h r o u g h experiment or t h e o r e t i c a l s t u d y . The t h i r d i s plasma s i m u l a t i o n t h r o u g h the use of computers.. By i t s v e r y n a t u r e , developments i n plasma s i m u l a t i o n have p a r a l l e l e d advances i n computer t e c h n o l o g y . From i t s crude b e g i n n i n g s a l i t t l e o ver twenty y e a r s ago, plasma s i m u l a t i o n has e v o l v e d i n t o a r i c h and d i v e r s e f i e l d of s t u d y . S i m u l a t i o n s a r e now a commonly used t o o l , complementing b o t h experiment and t h e o r y . I t can s e r v e as a way of i n d e p e n d e n t l y v e r i f y i n g r e s u l t s o b t a i n e d t h r o u g h o t h e r means. I t can a l s o be used t o examine problems t h a t might be i m p o s s i b l e t o a n a l y z e e x p e r i m e n t a l l y (due t o t e c h n o l o g i c a l c o n s t r a i n t s ) or a n a l y t i c a l l y (due t o m a t h e m a t i c a l c o m p l e x i t i e s ) . There have many d i f f e r e n t t y p e s of computer codes d e v e l o p e d . One of the more commonly used ones i s the P a r t i c l e I n C e l l (PIC) model. A program u t i l i z i n g t h e PIC scheme was used i n t h i s work. Another l a r g e f a m i l y of codes i n use a r e t h e Hydro codes. T h i s l a t t e r method models the hydrodynamic p r o c e s s e s o c c u r r i n g i n a plasma t h r o u g h the use of f l u i d e q u a t i o n s . 2 A q u e s t i o n the modeler must always ask h i m s e l f when a n a l y z i n g h i s computer g e n e r a t e d output i s : Does t h i s d ata r e f l e c t r e a l i t y ? W h i l e a program may run w i t h o u t g e n e r a t i n g e r r o r messages, and produce voluminous amounts of d a t a , t h i s i n i t s e l f does not ensure t h a t the d a t a i s not j u s t n u m e r i c a l g i b b e r i s h . T h i s i s e s p e c i a l l y a problem when the p h y s i c s b e i n g s i m u l a t e d has not been f u l l y e x p l o r e d e i t h e r by e x p e r i m e n t a l or t h e o r e t i c a l means. For i f t h i s i s the c a s e , t h e r e i s no independent means of v e r i f y i n g the v a l i d i t y of the program's r e s u l t s . The b e s t a l t e r n a t i v e then i s t o t e s t the program on known p h y s i c a l p r o c e s s e s and compare i t s r e s u l t s a g a i n s t the e s t a b l i s h e d t h e o r e t i c a l or e x p e r i m e n t a l r e s u l t s . I f the program a g r e e s w i t h the known r e s u l t s , i t i s a t l e a s t v a l i d f o r the t e s t e d c o n d i t i o n s and w i l l h o p e f u l l y s t i l l be p h y s i c a l l y v a l i d f o r the c o n d i t i o n s of i n t e r e s t ( i . e . new s i t u a t i o n s ) . T h i s i s the r a t i o n a l e of t h i s work. S i m u l a t i o n s of l a s e r plasma i n t e r a c t i o n s i n v o l v i n g the use of a PIC code have been p l a n n e d a t UBC, but b e f o r e t h i s work can be s t a r t e d however, the code t o be used must be t e s t e d . I t must be t e s t e d on a v a r i e t y of known phenomenon i n o r d e r t o see how w e l l i t behaves. I n t h i s work, an e l e c r o s t a t i c PIC code has been d e v e l o p e d and t e s t e d . R e s u l t s from s i m u l a t i o n s u s i n g t h i s code were compared t o the t h e o r y of Landau damping and the two stream i n s t a b i l i t y . 3 In c h a p t e r 2 t h e r e i s a b r i e f d i s c u s s i o n of the p e r t i n e n t plasma t h e o r y . T h i s i n c l u d e s a d e s c r i p t i o n of Debye s c r e e n i n g and plasma o s c i l l a t i o n s . A m a t h e m a t i c a l d e s c r i p t i o n of Landau damping i s n e x t , f o l l o w e d by a d i s c u s s i o n of the two stream i n s t a b i l i t y . C hapter 3 opens w i t h a an o v e r v i e w of the e v o l u t i o n of PIC codes. The second s e c t i o n d e s c r i b e s the plasma model s i m u l a t e d by the program. The nex t two s e c t i o n s o u t l i n e the scheme used by the program f o r moving the p a r t i c l e s . The c h a p t e r c o n c l u d e s w i t h a d i s c u s s i o n of t h r e e s c a l i n g p arameters and t h e i r s i g n i f i c a n c e . I n c h a p t e r 4 r e s u l t s of v a r i o u s s i m u l a t i o n s a r e p r e s e n t e d and d i s c u s s e d . In s e c t i o n 1 Landau damping i n plasma waves i s i n v e s t i g a t e d . S e c t i o n 2 d e a l s w i t h the two stream i n s t a b i l i t y . I n both s e c t i o n s the time e v o l u t i o n of energy ( k i n e t i c , e l e c t r o s t a t i c , and t o t a l ) i s examined. As w e l l the program i s t e s t e d f o r r e v e r s i b i l i t y . Comparisons w i t h the known t h e o r e t i c a l p r e d i c t i o n s a r e made. Good agreement i s o b t a i n e d , g i v i n g c o n f i d e n c e i n the r e l i a b i l i t y of the code. In c h a p t e r 5 some c o n c l u s i o n s a r e drawn. A d i s c u s s i o n of p o s s i b l e f u t u r e developments f o r t h i s code f o l l o w s . The two ap p e n d i c e s r e l a t e t o the code. In appendix A, the method of c a l c u l a t i n g the system energy i s g i v e n . Appendix B c o n s i s t s of documented l i s t i n g of the programs. 4 CHAPTER 2 THEORY I n t r o d u c t i o n In t h i s c h a p t e r o n l y a b r i e f r e v i e w w i l l be g i v e n of the plasma t h e o r y p e r t i n e n t t o t h i s work. T h i s w i l l a l s o s e r v e t o i n t r o d u c e c e r t a i n p h y s i c a l c o n s t a n t s of i n t e r e s t . Debye S c r e e n i n g The p o t e n t i a l produced by a p o i n t charge can be w r i t t e n as * ( r ) = q ( 2 . 1 ) r where q i s a t the o r i g i n . In a plasma however, because of charge m o b i l i t y , any bare charge w i l l be q u i c k l y shrouded by a c l o u d of o p p o s i t e c h a r g e s . T h i s w i l l t e n d t o n e u t r a l i z e any p o t e n t i a l t h a t might be seen by a n other charge o u t s i d e the s c r e e n i n g c l o u d . E q u i v a l e n t l y , the t e s t charge i n s i d e t h e c l o u d w i l l not see any p o t e n t i a l s from c h a r g e s beyond a c e r t a i n l e n g t h . Thus f o r d i s t a n c e s g r e a t e r than a s p e c i f i e d l i m i t , the p o t e n t i a l produced by a g i v e n charge i s assumed t o be n e g l i g i b l e . T h i s l i m i t i s c a l l e d the Debye l e n g t h and 5 i s d e f i n e d as (2.2) KT„ 4?rn e 2 o where K i s the Boltzmann's c o n s t a n t , T e the e l e c t r o n t e m p e r a t u r e , and n 0 the e l e c t r o n number d e n s i t y . For s e p a r a t i o n s g r e a t e r than the Debye l e n g t h , Coulomb p a r t i c l e i n t e r a c t i o n s can be n e g l e c t e d . T h i s e f f e c t i s known as Debye s c r e e n i n g . The Debye l e n g t h s e r v e s as s c a l e l e n g t h . Plasma phenomena o c c u r i n g over d i s t a n c e s g r e a t e r than X D (plasma waves f o r example), a r e c h a r a c t e r i z e d by c o l l e c t i v e p a r t i c l e b e h a v i o u r . In t h i s regime i t i s the motion of the c l o u d of p a r t i c l e s as a whole t h a t i s i m p o r t a n t . I n d i v i d u a l p a r t i c l e m o t i o n s do not p l a y a p a r t . However, f o r d i s t a n c e s of the o r d e r of X D or s h o r t e r , i n d i v i d u a l p a r t i c l e b e h a v i o u r d o m i n a t e s . Because of t h i s , waves w i t h wavelength £ X D become h e a v i l y damped. A c o n s t a n t t h a t w i l l be r e f e r r e d t o f r e q u e n t l y i n t h i s work i s the Debye wave number k Q. I t i s e q u a l t o X D " 1 . Plasma O s c i l l a t i o n s C o n s i d e r next what happens i n a c o l d plasma (KT=0), when t h e e l e c t r o n s a r e p e r t u r b e d as a group from t h e i r 6 e q u i l i b r i u m p o s i t i o n s . W i t h the i o n s f o r m i n g a f i x e d u n i f o r m background, the r e s u l t i n g charge s e p a r a t i o n s w i l l c r e a t e e l e c t r i c f i e l d s i n the plasma. These f i e l d s w i l l e x e r t a r e s t o r i n g f o r c e on the e l e c t r o n s . But because of th e i n e r t i a of t h e e l e c t r o n s t hey w i l l o s c i l l a t e back and f o r t h . The f r e q u e n c y of t h i s o s c i l l a t i o n i s g i v e n as HI, 47rn e 2 o (2.3) T h i s i s the plasma f r e q u e n c y . Plasma o s c i l l a t i o n s a r e seen as a c o l l e c t i v e phenomenon s i n c e i t i s the e l e c t r o n s as a whole t h a t o s c i l l a t e . A nother plasma c o n s t a n t t h a t s h o u l d be d e f i n e d i s t h e t h e r m a l v e l o c i t y . The t h e r m a l v e l o c i t y i s a measure of how hot the plasma i s and f o r one d i m e n s i o n i t i s d e f i n e d as KT e m (2.4) E q u a t i o n s ( 2 . 2 ) , (2.3) and (2.4) can be combined t o g i v e the r e l a t i o n = v. (2.5) I f the plasma has a non-zero temperature the t h e r m a l motions of the e l e c t r o n s a l l o w f o r the t r a n s f e r of 7 i n f o r m a t i o n from one r e g i o n of the plasma t o a n o t h e r . As a r e s u l t , CJ becomes a f u n c t i o n of the wave number k. The group v e l o c i t y i s no l o n g e r z e r o and so p e r t u r b i n g the e l e c t r o n s w i l l r e s u l t i n plasma waves b e i n g g e n e r a t e d . The r e l a t i o n s h i p between w and k b r i n g s us t o the t o p i c of the next s e c t i o n , Landau damping. Landau Damping As Landau damping p l a y s an i n t e g r a l r o l e i n the p h y s i c s s i m u l a t e d i n t h i s work, some mention of the u n d e r l y i n g t h e o r y , f i r s t d e s c i b e d by Landau(1946), s h o u l d be made. T h e r e f o r e i n t h i s s e c t i o n , a b r i e f t h e o r e t i c a l d e s c r i p t i o n of l i n e a r Landau damping i s g i v e n . The d i s c u s s i o n w i l l be l i m i t e d t o the one d i m e n s i o n a l problem t a c k l e d i n the s i m u l a t i o n s . In o r d e r t o p r o c e e d , two e q u a t i o n s a r e r e q u i r e d . They a r e P o i s s o n ' s e q u a t i o n = -4irp (2.6) or a l t e r n a t i v e l y 8 V.E m -4ire y f ( v ) d v (2.7) and Boltzmann's e q u a t i o n . 3 f + w f + F 3 f 3 t m 3v 3f 3t (2.8) Where = e E + ^ x B (2.9) Assume t h a t we have a q u a s i - n e u t r a l homogeneous plasma. F u r t h e r m o r e , assume t h a t the e l e c t r o n s are n o n - r e l a t i v i s t i c ( v t h « c ) so t h a t F = qE, and t h a t t h e i o n s a r e f i x e d . Thus, the i o n d i s t r i b u t i o n f u n c t i o n remains c o n s t a n t and i s e q u a l t o n 0 , t h e plasma d e n s i t y . The plasma w i l l a l s o be assumed t o be c o l l i s i o n l e s s w i t h i n the t i m e frame of i n t e r e s t . T h e r e f o r e , the c o l l i s i o n term i n e q u a t i o n (2.8) O f / 3 t ) | =0 The e l e c t r o n d i s t r i b u t i o n i s p e r t u r b e d a t t=0 so t h a t f ( x , v , 0 ) = f 0 ( v ) + f , ( x , v , 0 ) F o r l i n e a r damping, f 0 » f , , where f 0 ( v ) i s t h e u n p e r t u r b e d p a r t of t h e d i s t r i b u t i o n . I f t h e r e a r e no e x t e r n a l l y a p p l i e d f i e l d s , as was the case h e r e , any f i e l d s p r e s e n t i n t h e plasma a r e due t o f l u c u a t i o n s i n the e l e c t r o n d e n s i t y r e p r e s e n t e d by f , . E i s t h u s a f i r s t o r d e r term. W i t h t h e s e assumptions i n mind, e q u a t i o n (2.8) t o f i r s t o r d e r , 9 can be w r i t t e n as 3f af —- 1 + v l - e E 3 f 3 t 3x m 1 g (2.10) E q u a t i o n (2.10) i s the l i n e a r i z e d V l a s o v e q u a t i o n . I f the p e r t u r b a t i o n c o n s i s t s of a one d i m e n s i o n a l p l a n e wave t r a v e l l i n g i n the x d i r e c t i o n , f , and E 1 t a k e the form A • exp( i ( k x - o t ) ) . The wave f r e q u e n c y CJ can be e i t h e r r e a l or complex. E q u a t i o n (2.10) now becomes - i u f + i k f v • e E 3 f 3v or f = i e E 9 f 0 1 1 m I - (2.11) 3v kv - u Combining e q u a t i o n s (2.7) and (2.11) we get , u>2 r 3f i 1 = - I I o d v kn J 3v k v - u (2.12) w i t h w p, t h e plasma f r e q u e n c y , d e f i n e d by e q u a t i o n (2.3) For a M a x w e l l i a n plasma f 0 = n 0 f ( v ) where f ( v ) i s the n o r m a l i z e d M a x w e l l i a n . E q u a t i o n 2.10 can then be w r i t t e n as 10 1 = -a) n 2 f 3f. 1 dv k J 3v k v - u (2.13) Equat i o n (2.13) i s the d i s p e r s i o n r e l a t i o n . A s t r a i g h t - f o r w a r d e v a l u a t i o n of the i n t e g r a l i s not p o s s i b l e due t o the p o l e a t v=cj/k. N o r m a l l y , i n t e g r a l s of t h i s type a r e s o l v e d u s i n g the r e s i d u e theorem so t h a t , where C, i s a c o n t o u r a l o n g the r e a l a x i s c o n t o u r C 2 i s a s e m i c i r c l e of i n f i n i t e r a d i u s i n one of the i m a g i n a r y h a l f p l a n e s and R(w/k) i s the r e s i d u e a t the p o l e v=w/k. For t h i s method t o work the i n t e g r a l f o r c o n t o u r C 2 must v a n i s h f o r v—>±i». U n f o r t u n a t e l y f o r the case where W(v)= ( 9 f 0 / 9 v ) [ v - w / k ] " 1 , the C 2 i n t e g r a l blows up i n e i t h e r h a l f p l a n e as 9 f 0 / 3 v c o n t a i n s the term exp(-v2/v.„ 2 ) . F o r t h t h i s reason the i n t e g r a l must e v a l u a t e d n u m e r i c a l l y . C a l c u l a t e d t h e o r e t i c a l v a l u e s f o r u/u and f o r s e v e r a l v a l u e s of k/k D a r e g i v e n i n T a b l e I , where 7 i s the im a g i n a r y p a r t of the wave f r e q u e n c y . These numbers were taken from C a n o s a ( 7 3 ) . (2.14) 11 T a b l e I T h e o r e t i c a l V a l u e s For Landau Damping And F r e q u e n c i e s -7A> p 0.2 1 .0640 5.511E-05 0.3 1.1598 1.262E-02 0.4 1.2850 0.0661 0.5 1 .4146 0.1534 0.6 1.5457 0.2641 0.8 1.7999 0.5346 1 .0 2.0459 0.8513 2.0 3.1891 2.8272 The Two Stream I n s t a b i l i t y A nother e l e c t r o s t a t i c phenomenom t h a t was s i m u l a t e d i s the two st r e a m i n s t a b i l i t y . I n t h i s case t h e u n p e r t u r b e d plasma c o n s i s t s of two e l e c t r o n streams of e q u a l d e n s i t y moving t h r o u g h a f i x e d i o n background w i t h e q u a l and o p p o s i t e v e l o c i t i e s . F i s then w r i t t e n as 12 f (v) = n 0 «(v - v ) + n 0 6(V + V ) (2.15) where v 0 i s the s t r e a m i n g v e l o c i t y . The plasma i s then p e r t u r b e d i n a manner s i m i l a r t o t h e p r o c e d u r e of the p r e v i o u s s e c t i o n . To d e s c r i b e t h i s p r o c e s s m a t h e m a t i c a l l y , we can p r o c e e d d i r e c t l y from e q u a t i o n ( 2 . 1 2 ) , so t h a t , 1 = _P 2 [* , ( V- V dv + f J kv - u J *'<v + V dv kv - to To s o l v e the i n t e g r a l s we make use of t h e i d e n t i t y (2.16) 6* (x - a)f (x)dx • - f • (a) (2.17) By making the a p p r o p r i a t e s u b s t i t u t i o n s , we get CO' 1 = _P_ 2 (kv - oj) 2 0 (kv + w)2 0 (2.18) T h i s i s the d i s p e r s i o n r e l a t i o n f o r the two stream i n s t a b i l i t y . E q u a t i o n (2.18) can be r e a r r a n g e d so t h a t , - w 2 [ 2 ( k v 0 ) 2 + w p 2 ] + ( k v 0 ) 2 [ ( k v 0 ) 2 - w D 2 ] = 0 (2.19) There w i l l be f o u r r o o t s t o the e q u a t i o n . The u2 r o o t s a r e g i v e n by 13 U 2 . ( K , 2 + ^ ± ^ [ 8(kV )» + U P 2 ] H ( 2 > 2 0 ) I f kv 0>cj p then <j2 w i l l a l ways be p o s i t i v e and a l l f o u r r o o t s w i l l be r e a l . There w i l l be no. i n s t a b i l i t y and the d i s t r i b u t i o n f u n c t i o n w i l l remain b a s i c a l l y unchanged. I f kv 0< oip t h e r e w i l l be two r e a l r o o t s and two im a g i n a r y c o n j u g a t e r o o t s s i n c e 2 * u 2 = ( k v , 2 + • [ B O O M 2 + U , 2 ] ° 2 2 < 0 f o r k v 0 < cop (2.21 ) One of the im a g i n a r y r o o t s w i l l r e s u l t i n a growing i n s t a b i l i t y s i n c e f w i l l then c o n t a i n an e x p o n e n t i a l term t h a t grows w i t h t i m e . I f ± 7 are the i m a g i n a r y r o o t s , then t h e maximum v a l u e f o r \y\ o c c u r s k v 0 / c j = /3/8 and i s e q u a l t o I7 | =u //8~~ . " max P V a l u e s f o r | y/cop \ are g i v e n i n T a b l e I I . 14 T a b l e I I T h e o r e t i c a l Growth Rates And F r e q u e n c i e s For The Two Stream I n s t a b i l i t y kv 0/w p w/cj p \t\/< 0.0 1 .0000 0.0000 0.1 1 .0147 0.0981 0.2 1 .0557 0.1856 0.3 1.1161 0.2564 0.4 1 . 1895 0.3082 0.5 1.2712 0.3406 0.6 1 .3583 0.3534 0.7 1.4488 0.3450 0.8 1.5417 0.3113 0.9 1.6363 0.2398 1.0 1.7320 0.0000 15 CHAPTER 3 DESCIPTION OF PIC CODES I n t r o d u c t i o n For a r e v i e w of PIC codes i n g e n e r a l , t h e reader i s r e f e r r e d t o the r e v i e w a r t i c l e s by Dawson (83) or Langdon (76) or the book by Hockney ( 8 1 ) . At i t s most b a s i c l e v e l a p a r t i c l e code keeps t r a c k of a c o l l e c t i o n of charged p a r t i c l e s t h a t move about a c c o r d i n g t o t h e i r s e l f - c o n s i s t e n t e l e c t r i c f i e l d s . At f i r s t g l a n c e t h i s might seem t o be the i d e a l s o l u t i o n , r e c r e a t e the plasma environment down t o the l a s t e l e c t r o n . Of c o u r s e t h i s approach i s h i g h l y u n f e a s i b l e . As an i l l u s t r a t i o n of how u n f e a s i b l e i t i s , assume t h a t we wanted t o s i m u l a t e a plasma t h a t c o n t a i n e d 10 1° p a r t i c l e s . T h i s i s a s m a l l number when compared t o t h e number of p a r t i c l e s found i n l a b o r a t o r y or n a t u r a l plasmas. F u r t h e r m o r e assume we have a computer t h a t can p e r f o r m an o p e r a t i o n i n 1 0 " 1 0 seconds, a speed w e l l beyond p r e s e n t day l i m i t s . To c a l c u l a t e t h e a c c e l e r a t i o n of each p a r t i c l e i n a plasma of n p a r t i c l e s a t l e a s t n 2 c a l c u l a t i o n s must be made. So, f o r 10 1° p a r t i c l e s i t would t a k e a minimum of 300 y e a r s t o s i m u l a t e one time s t e p . A normal run might need anywhere up t o a few thousand time s t e p s . 16 As any r e a l plasma w i l l have a p a r t i c l e number t h a t i s many o r d e r s of magnitude g r e a t e r than' can p o s s i b l y be f o l l o w e d i n a p r a c t i c a l s i m u l a t i o n , a compromise must be made. The normal p r a c t i c e i s t o make one s i m u l a t e d p a r t i c l e e q u i v a l e n t t o a group of p a r t i c l e s i n a r e a l plasma. So, i n moving one s i m u l a t e d p a r t i c l e , one might a c t u a l l y be moving hundreds or thousands of e l e c t r o n s or i o n s a l l w i t h the same v e l o c i t y . There s t i l l remains the problem of c a l c u l a t i n g the f o r c e a c t i n g on each p a r t i c l e . E a r l y codes d i d t h i s by summing up a l l the f o r c e s a c t i n g on a p a r t i c l e (see Dawson(62)). In a plasma of n p a r t i c l e s t h i s would r e q u i r e of t h e o r d e r of n 2 c a l c u l a t i o n s . T h i s n 2 f a c t o r l i m i t s t h e s i z e of s i m u l a t i o n s , so t h a t runs w i t h many thousands of p a r t i c l e s a r e not p r a c t i c a l u s i n g t h i s method. For t h i s r e a s o n t h e c o n c e p t s of a mesh and charge s h a r i n g were d e v e l o p e d . A f i x e d s e t of s p a t i a l g r i d p o i n t s , or mesh, i s p l a c e d i n the plasma. A s s i g n e d t o t h e s e g r i d p o i n t s a r e v a l u e s f o r the Coulomb e l e c t r i c f i e l d and the p o t e n t i a l d e t e r m i n e d by the p a r t i c l e p o s i t i o n s w i t h i n the g r i d . The f o r c e a c t i n g on a p a r t i c l e i s then found from the v a l u e s on the mesh p o i n t s . The f i r s t such scheme was c a l l e d the N e a r e s t G r i d P o i n t (NGP) method. I t was d e v e l o p e d a t S t a n f o r d i n d e p e n d e n t l y f o r 1-D by Burger (65) and f o r 2-D by Hockney ( 6 6 ) . In t h i s 17 scheme, when c a l c u l a t i n g the p o t e n t i a l , the g r i d p o i n t a p a r t i c l e was n e a r e s t t o , r e c e i v e d a l l i t s c h a r g e . Though c r u d e , i t d i d p r o v i d e u s e f u l r e s u l t s and i t was c o n s i d e r a b l y f a s t e r than p r e v i o u s methods. A f u r t h e r r e f i n e m e n t b r i n g s us t o the type of code used i n t h i s work. I n s t e a d of a s s i g n i n g the t o t a l p a r t i c l e charge t o j u s t one g r i d p o i n t , the charge i s sh a r e d by the two mesh p o i n t s between which the p a r t i c l e l i e s . Charge s h a r i n g was f i r s t d e v e l o p e d by B i r d s a l l and Fuss ( 6 9 ) . The p a r t i c l e c o u l d now be seen as b e i n g s p r e a d or d i f f u s e d over the w i d t h of two c e l l s (a c e l l b e i n g the d i s t a n c e between two mesh p o i n t s ) , a n a l o g o u s t o the i d e a of a c h a r g e , c l o u d . Because of t h i s , t h i s scheme was c h r i s t e n e d the C l o u d In C e l l (CIC) method. The c o n t r i b u t i o n t o t h e charge a t b o t h mesh p o i n t s i s de t e r m i n e d by l i n e a r i n t e r p o l a t i o n . The charge d e n s i t y a t each mesh p o i n t i s found by summing over a l l the p a r t i c l e s . The e l e c t r i c f i e l d a t t h e mesh p o i n t s i s then c a l c u l a t e d from P o i s s o n ' s e q u a t i o n . I n o r d e r t o o b t a i n the new v e l o c i t y of a p a r t i c l e , t he e l e c t r i c f i e l d E , a t the p o s i t i o n of the p a r t i c l e i s r e q u i r e d . To a c h i e v e t h i s the same i n t e r p o l a t i o n scheme i s used, but i n a r e v e r s e f a s h i o n . The e l e c t r i c f i e l d a t b o t h mesh p o i n t s on e i t h e r s i d e of the p a r t i c l e c o n t r i b u t e s t o the v a l u e of the e l e c t r i c f i e l d a t the p a r t i c l e . More a c c u r a t e schemes have been d e v e l o p e d where the 18 p a r t i c l e charge i s spread over more than one c e l l . However a p r i c e i s p a i d f o r the i n c r e a s i n g a c c u r a c y i n the form of i n c r e a s i n g c o m p u t a t i o n t i m e . For a d i s c u s s i o n of t h e s e h i g h e r o r d e r methods, the r e a d e r i s r e f e r r e d t o Hockney ( 8 1 ) . A common g e n e r i c name f o r t h e s e t y p e s of schemes i s , P a r t i c l e i n C e l l (PIC) code. PIC codes have been de v e l o p e d t o s t u d y many d i f f e r e n t s i t u a t i o n s . For l a s e r - p l a s m a i n t e r a c t i o n s and r e l a t i v i s t i c plasmas, t h e r e are e l e c t r o - m a g n e t i c (EM) PIC codes. These programs, as w e l l as c a l c u l a t i n g t h e Coulomb f i e l d , d e t e r m i n e the e l e c t r o - m a g n e t i c E and B f i e l d s on the mesh p o i n t s . For c a s e s where t h e r e a r e no a p p l i e d f i e l d s , t h e r e a r e the E l e c t r o s t a t i c (ES) PIC codes where o n l y the Coulomb f i e l d s a r e c a l c u l a t e d . PIC codes have been w r i t t e n f o r 1, 2 , and 3 d i m e n s i o n a l s i m u l a t i o n s . The methods used do not r e a l l y change t h a t much w i t h a change i n d i m e n s i o n . E a r l y e f f o r t s c o n s i s t e d of s o l e l y 1 d i m e n s i o n a l s i m u l a t i o n s as t h e s e were the o n l y t y p e s f e a s i b l e g i v e n the l e v e l of computer t e c h n o l o g y t h a t e x i s t e d a t t h a t t i m e . S i n c e t h e n , the t e c h n o l o g y has improved g r e a t l y , but d e s p i t e t h i s , 3 d i m e n s i o n a l s i m u l a t i o n s a r e s t i l l r a r e . S i m u l a t i o n s i n 2 d i m e n s i o n s a r e more common, but even today, much of the plasma s i m u l a t i o n s a r e done i n one d i m e n s i o n . The problem l i e s i n t h a t the number of p a r t i c l e s r e q u i r e d t o do e q u i v a l e n t s i m u l a t i o n s grows g e o m e t r i c a l l y w i t h the number of d i m e n s i o n s used. As 19 w e l l , the number of q u a n t i t i e s t o c a l c u l a t e grows s u b s t a n t i a l l y w i t h an i n c r e a s e i n d i m e n s i o n . So, even w i t h the most advanced of computers, 3 d i m e n s i o n a l s i m u l a t i o n s have remained r a t h e r c o a r s e l y g r a i n e d . There has been more s u c c e s s w i t h 2 d i m e n s i o n a l codes. They have been u s e f u l i n a n a l y s i n g many p r o c e s s e s t h a t o c c u r i n plasmas p a r t i c u l a r l y t h o s e such as two plasmon decay t h a t cannot be s i m u l a t e d i n a 1 d i m e n s i o n a l model because of t h e i r geometry. However t h e r e i s s t i l l a g r e a t d e a l l e f t t o e x p l o r e t h a t can be h a n d l e d i n a 1 d i m e n s i o n a l model. The program used i n t h i s work was a 1 d i m e n s i o n a l PIC code. The Model The type of model t h i s program s i m u l a t e s i s one of the s i m p l e s t . The plasma m o d e l l e d i s an i n f i n i t e one d i m e n s i o n a l plasma, w i t h p e r i o d i c boundary c o n d i t i o n s . I t i s homogeneous w i t h a f i x e d i o n background and t h e r e a r e no a p p l i e d f i e l d s . I t i s a l s o a c o l l i s i o n l e s s model. What f o l l o w s i s a more d e t a i l e d d e s c r i p t i o n of t h e s e v a r i o u s f e a t u r e s . The most o b v i o u s f e a t u r e of a one d i m e n s i o n a l model i s t h a t q u a n t i t i e s can v a r y a l o n g o n l y one a x i s . To be c o n s i s t e n t w i t h the n o t a t i o n used i n t h i s work, l e t t h i s a x i s be the x a x i s . Charge d e n s i t y t h e n , v a r i e s o n l y i n the 20 x d i r e c t i o n and i s assumed c o n s t a n t i n the y and z d i r e c t i o n . C o n s e q u e n t l y the e l e c t r i c f i e l d g e n e r a t e d from the v a r i a t i o n s i n the charge s e p a r a t i o n v a r i e s o n l y i n the x d i r e c t i o n . From t h i s i t f o l l o w s t h a t p a r t i c l e s w i l l undergo a c c e l e r a t i o n s o n l y i n the x d i r e c t i o n . P a r t i c l e v e l o c i t i e s i n the y and z d i r e c t i o n s w i l l be c o n s t a n t . There i s something here t h a t s h o u l d be c l a r i f i e d . The charge d e n s i t y and t h e n e u t r a l plasma d e n s i t y a r e two d i s t i n c t q u a n t i t i e s . The charge d e n s i t y i s a measure of the e x c e s s e l e c t r o n charge i n a c e r t a i n volume. For r e g i o n s where the e l e c t r o n charge exceeds the i o n i c charge t h e r e e x i s t s a n e g a t i v e charge d e n s i t y . S i m i l a r l y f o r r e g i o n s where the o p p o s i t e c o n d i t i o n s p r e v a i l , t h e r e i s a p o s i t i v e charge d e n s i t y . For n e u t r a l r e g i o n s t h e r e i s a z e r o charge d e n s i t y . I f we assume t h a t the plasma i s q u a s i - n e u t r a l , meaning, i f we n e g l e c t s m a l l p e r t u r b a t i o n s , n =n =n 0, then n 0 i s c o n s i d e r e d t o be the n e u t r a l plasma d e n s i t y . N and n a r e r e s p e c t i v e l y , the e l e c t r o n and i o n number d e n s i t i e s . S i n c e the plasma i s i n f i n i t e , we d i d not have t o use boundary c o n d i t i o n s t h a t would r e f l e c t the p r o c e s s e s o c u r r i n g a t the a c t u a l or p h y s i c a l edge of a plasma. These type of boundary c o n d i t i o n s g e n e r a l l y i n t r o d u c e c o m p l i c a t i o n s t o t h e program. W i t h an i n f i n i t e plasma however, t h i n g s a r e s i m p l e r . To i n c o r p o r a t e p e r i o d i c boundary c o n d i t i o n s you o n l y have t o ensure t h a t the charge d e n s i t i e s a t b o t h b o u n d a r i e s a r e e q u a l t o each o t h e r . There 21 s t i l l remains the problem of what t o do w i t h p a r t i c l e s t h a t l e a v e the plasma r e g i o n . The method we used i s q u i t e s t a n d a r d , Langdon(76). When a p a r t i c l e e x i t s from one s i d e w i t h a v e l o c i t y v i t i s r e - i n t r o d u c e d i n t o the plasma a t the o t h e r s i d e w i t h the same v e l o c i t y v. A v a r i a t i o n on t h i s method i s t o a s s i g n the p a r t i c l e a new random v e l o c i t y . T h i s v e l o c i t y i s g e n e r a t e d from the p r o p e r h a l f of a M a x w e l l i a n d i s t r i b u t i o n , so as t o keep the p a r t i c l e moving i n the same d i r e c t i o n as b e f o r e . The prime m o t i v e f o r u s i n g t h i s method i s t h a t i t r e f l e c t s the p e r i o d i c n a t u r e of the b o u n d a r i e s . I t i s a l s o an easy way of e n s u r i n g t h a t the n e u t r a l d e n s i t y remains c o n s t a n t . T h i s b r i n g s us on t o the next p o i n t . As was s t a t e d e a r l i e r , the plasma m o d e l l e d i s homogeneous. The n e u t r a l d e n s i t y i s c o n s t a n t t hroughout the plasma r e g i o n . T h i s i s i n c o n t r a s t t o some EM PIC codes where a d e n s i t y ramp i s used. T h i s model a l s o d i f f e r s from an EM code i n t h a t t h e r e a r e no a p p l i e d f i e l d s (There i s no l a s e r s h i n i n g on the p l a s m a ) . The o n l y e l e c t r i c f i e l d s p r e s e n t a r e t h o s e g e n e r a t e d by charge s e p a r a t i o n s . T h i s i s why t h i s type of code i s g e n e r a l l y known as an e l e c t r o s t a t i c model. T h i s does not i m p l y t h a t the c h a r g e s a r e f i x e d . The model i s o n l y c a l l e d t h i s so as t o d i f f e r e n t i a t e i t from models t h a t do i n c o r p o r a t e e l e c t r o - m a g n e t i c f i e l d s . The plasma i s assumed t o be n o n - r e l a t i v i s t i c (V t h /C«1 ) so t h a t t h e r e a r e 22 n e g l i g i b l e magnetic f o r c e s ( [ v / c ] x B) p r e s e n t . In t h e program o n l y t h e e l e c t r o n s a r e moved. The i o n s a r e immobile, so they a r e seen as c o n s t a n t p o s i t i v e background c h a r g e . There i s no problem w i t h t h i s as l o n g as the time s c a l e s a r e much s h o r t e r than u ' 1 , where o> i s the i i i o n f r e q u e n c y As time d u r a t i o n s f o r runs i n t h i s work were of the o r d e r 50a/ 1 which i s much l e s s than w " 1 , no d i f f i c u l t i e s a r o s e from k e e p i n g the i o n s f i x e d . F i n a l l y we note t h a t t h i s i s a c o l l i s i o n l e s s plasma. The model does not ta k e i n t o account any e f f e c t s t h a t would be due t o p a r t i c l e c o l l i s i o n s . T h i s assumption i s v a l i d as l o n g as the e l e c t r o n c o l l i s i o n f r e q u e n c y i s much lower than the e l e c t r o n o s c i l l a t i o n f r e q u e n c y . S i n c e v c 1 ln(12Trn X 3) a 0 a> p as l o n g as TIX* » 1, c o l l i s i o n a l e f f e c t s can be i g n o r e d . 23 N o t a t i o n B e f o r e p r o c e e d i n g w i t h the d i s c u s s i o n of the p a r t i c l e moving scheme, a mention s h o u l d be made of the n o t a t i o n used. The n o t a t i o n used i n the v a r i o u s d i f f e r e n c i n g e q u a t i o n s i s f a i r l y c o n v e n t i o n a l . The v a r i o u s q u a n t i t i e s c a l c u l a t e d , such as the e l e c t r i c f i e l d , p a r t i c l e p o s i t i o n and so on, a r e marked by s u p e r s c r i p t and s u b s c r i p t c h a r a c t e r s . The s u p e r s c r i p t s n, n+1/2, and n+1 f o r any q u a n t i t y denote t h e ti m e f o r wh i c h t h a t v a l u e i s d e t e r m i n e d . When we t a l k of time n, the time r e f e r r e d t o i s a c t u a l l y e q u a l t o n«At, where n i s the number of time s t e p s and At i s the time increment per time s t e p . The s u b s c r i p t s f o r the charge d e n s i t y and the e l e c t r i c f i e l d denote t h e mesh p o i n t a t which the q u a n t i t y i s d e t e r m i n e d . These s u b s c r i p t s can be e i t h e r the l e t t e r s i or j . The p a r t i c l e p o s i t i o n and v e l o c i t y have a s u b s c r i p t p. T h i s s e r v e s as a reminder t h a t p o s i t i o n and v e l o c i t y a r e r e l a t e d d i r e c t l y t o the p a r t i c l e and a r e not a s s o c i a t e d w i t h any s p a t i a l mesh. Other forms of s u p e r s c r i p t and s u b s c r i p t n o t a t i o n found i n t h i s work f o l l o w s i m i l a r p a t t e r n s . At t h i s p o i n t t h e r e a r e t h r e e c o n s t a n t s t h a t s h o u l d be i n t r o d u c e d . The number of mesh p o i n t s i n the s p a t i a l g r i d w i l l be g i v e n by N M . The t o t a l number of p a r t i c l e s i n a s i m u l a t i o n w i l l be denoted by N p . The number of p a r t i c l e s per c e l l w i l l be w r i t t e n as N c ( = N p / [ N M - 1 ] ) . 24 The Scheme In t h i s s e c t i o n a d e s c r i p t i o n of the program proc e d u r e f o r one time s t e p w i l l be g i v e n . The g e n e r a l scheme i s s t r a i g h t f o r w a r d and i s i l l u s t r a t e d i n f l o w c h a r t form i n F i g u r e 1. At the s t a r t of t i m e s t e p n the e l e c t r i c f i e l d s f o r time n a t a l l mesh p o i n t s a r e known. As w e l l , the p o s i t i o n s of the p a r t i c l e s a t time n and the v e l o c i t i e s of t h e p a r t i c l e s a t time n-1/2 a r e known. The f i r s t s t e p i s t o c a l c u l a t e the e l e c t r i c f i e l d a t the i n d i v i d u a l p a r t i c l e p o s i t i o n . T h i s i s done by f i r s t d e t e r m i n i n g between which two mesh p o i n t s the p a r t i c l e l i e s . U s i n g the e l e c t r i c f i e l d v a l u e s a t the s e two mesh p o i n t s , a v a l u e f o r E a t the p a r t i c l e p o s i t i o n ( E p ) i s found by l i n e a r i n t e r p o l a t i o n . The degree t o which a mesh p o i n t c o n t r i b u t e s t o E i s d e t e r m i n e d by a w e i g h t i n g term. T h i s term v a r i e s l i n e a r l y from a maximum of 1 when the p a r t i c l e i s e x a c t l y on the mesh p o i n t t o a minimum of 0 when the p a r t i c l e i s a mesh l e n g t h away. I f a p a r t i c l e i s between mesh p o i n t s j and j+1, E can be w r i t t e n as n n n E = a • E + • [ 1 - a] • E 0 £ a £ 1 p j j + 1 (3.1) As t h e mesh l e n g t h i s of t h e o r d e r of t h e Debye wavelength, Scheme For P a r t i c l e Mover F i g u r e 1 26 we do not c o n s i d e r e l e c t r i c f i e l d c o n t r i b u t i o n s from a d d i t i o n a l mesh p o i n t s . These c o n t r i b u t i o n s are deemed t o be n e g l i b l e due t o Debye s c r e e n i n g . Once E p has been c a l c u l a t e d we can now f i n d the p a r t i c l e v e l o c i t y f o r time n+1/2. S i n c e the f o r c e a c t i n g on the p a r t i c l e i s F = m'A = q' • E , (3.2) and the a c c e l e r a t i o n can be e x p r e s s e d i n d i f f e r e n c e form as n+V2 n-V 2 V + V p p n A f = A (3.3) the p a r t i c l e v e l o c i t y a t time n+1/2 can be w r i t t e n as n+ 1/ 2 n - 1 / 2 n V = V + g'.At-E p p m' p (3.4) Note t h a t i n d e t e r m i n i n g V p n *"2 we used the c e n t r a l d i f f e r e n c i n g scheme. The advantage of t h i s method over e i t h e r f o r w a r d or backward d i f f e r e n c i n g i s t h a t the former i s a c c u r a t e t o second o r d e r w h i l e the l a t t e r two a r e o n l y a c c u r a t e t o f i r s t o r d e r . We mark th e p a r t i c l e mass term m and t h e p a r t i c l e charge term q w i t h a p o s t r o p h e s as a reminder t h a t t h e s e model p a r t i c l e s r e p r e s e n t s h e e t s of e l e c t r o n s r a t h e r than a 27 s i n g l e e l e c t r o n . In e q u a t i o n (3.4) though, q'/m' can be r e p l a c e d by e/m, the e l e c t r o n charge t o mass r a t i o , s i n c e q' and e a r e r e l a t e d t o each o t h e r by the same s c a l e f a c t o r as m' and m a r e . As shown l a t e r , w h e n c a l c u l a t i n g some of the d i a g n o s t i c q u a n t i t i e s t h i s s c a l e f a c t o r does become i m p o r t a n t . Once we have V a t time n+1/2, the p a r t i c l e p o s i t i o n a t time n+1,X n + 1 , can be found. U s i n g c e n t r a l d i f f e r e n c e s a g a i n we get n+1 n n+V2 X = X + At«V P P P (3.5) Once a l l t h e new p o s i t i o n s have been d e t e r m i n e d , the charge d e n s i t y p f o r time n+1 can be found a t the mesh p o i n t s . To do t h i s , the same w e i g h t i n g scheme t h a t was used i n c a l c u l a t i n g the e l e c t r i c f i e l d a t the p a r t i c l e p o s i t i o n i s used. I f the p a r t i c l e i s between mesh p o i n t s j and j+1 a p o r t i o n of the p a r t i c l e charge i s g i v e n t o p " + 1 and the r e s t goes t o P " * ] n+1 n+1 P = P + a«qdx (3.6) j j n+1 n+1 P = P + [ l - a ] - q d x (3.7) j+1 j+1 28 where the v a l u e f o r a i s de t e r m i n e d i n the same way i t was when f i n d i n g E p . The c o n s t a n t qdx c o r r e s p o n d s t o a s i n g l e p a r t i c l e charge t a k i n g i n t o account the d e s i r e d t r u e plasma d e n s i t y . I t i s e q u a l t o qdx = e-n 0/N (3.8) c As the model i s supposed t o r e p r e s e n t an i n f i n i t e one d i m e n s i o n a l plasma, s p e c i a l c o n s i d e r a t i o n must be g i v e n t o the charge d e n s i t y a t the boundary mesh p o i n t s . I n t h i s c a s e , s i n c e p e r i o d i c boundary c o n d i t i o n s were used, the charge d e n s i t y a t e i t h e r boundary must be e q u a l . To a c h i e v e t h i s , a c o n t r i b u t i o n t o p a t e i t h e r boundary i s added t o a r u n n i n g t o t a l . A f t e r i t e r a t i n g t h r o u g h a l l t h e p a r t i c l e s t h i s t o t a l i s a s s i g n e d t o p a t b o t h b o u n d a r i e s . Once t h i s i s done, the c o n s t a n t i o n charge d e n s i t y must be s u b t r a c t e d from p t o g i v e us the net charge d e n s i t y (-Nc«qdx) . The next s t e p i s t o o b t a i n the v a l u e f o r p h a l f w a y between mesh p o i n t s . To do t h i s we t a k e the s i m p l e average n+1 n+1 n+1 P = [p +' P 1/2 J + 1 / 2 j j + 1 (3.9) We do t h i s so t h a t we can use c e n t r a l d i f f e r e n c e s t o o b t a i n t h e e l e c t r i c f i e l d a t t i m e n+1. The e l e c t r i c f i e l d i s c a l c u l a t e d from P o i s s o n ' s e q u a t i o n . In d i f f e r e n c e form t h i s can be w r i t t e n as 29 n+1 n+1 n E E + 4 77 • Ax • p j+1 (3.10) To c o r r e c t f o r any o f f s e t t h a t might o c c u r f o r n u m e r i c a l r e a s o n s , or o t h e r w i s e , a c o r r e c t i o n term i s s u b t r a c t e d from the e l e c t r i c f i e l d a t each mesh p o i n t . T h i s term i s found by summing up the e l e c t r i c f i e l d over a l l the mesh p o i n t s and d i v i d i n g by the number of mesh p o i n t s . I f t h e r e i s no o f f s e t t h i s term s h o u l d be z e r o . At t h i s p o i n t we have the e l e c t r i c f i e l d and p a r t i c l e p o s i t i o n s f o r time n+1 and the p a r t i c l e v e l o c i t i e s a t time n+1/2. The c y c l e i s now com p l e t e . I f t h i s i s not the f i n a l t ime s t e p , the c y c l e r e p e a t s , a f t e r f i r s t b l a n k i n g out the p a r r a y . O t h e r w i s e the program s w i t c h e s o u t , does some f i n a l d i a g n o s t i c s , i f any, and s t o p s . As a f i n a l n o t e , t h e r e s h o u l d be a b r i e f d e s c r i p t i o n g i v e n of the d a t a p r e p a r a t i o n t h a t i s n e c c e s s a r y b e f o r e the d i a g n o s t i c s can be performed. L e t ' s say the program was t o a n a l y s e the d a t a a t the end of a time s t e p u s i n g as i n p u t t h e v e l o c i t i e s and the e l e c t r i c f i e l d as c a l c u l a t e d above. The e l e c t r i c f i e l d v a l u e s would be those f o r time n+1, w h i l e t h e v e l o c i t i e s would be the v e l o c i t i e s f o r time n+1/2. For some d i a g n o s t i c r o u t i n e t h a t uses or compares both t h e s e q u a n t i t i e s t h e n , an e r r o r might be i n t r o d u c e d . For v e r y s h o r t time s t e p s , t h i s might not be a problem but f o r l a r g e r t i m e s t e p s i t w i l l become more of a f a c t o r . To a v o i d t h i s 30 problem, we made sure a l l the d a t a used by the d i a g n o s t i c package was d a t a f o r time n. T h i s e n t a i l e d p l a c i n g a t time n, the e l e c t r i c f i e l d , the p a r t i c l e p o s i t i o n s and v e l o c i t i e s r e q u i r e d , and the non-mesh-centered charge d e n s i t y i n s t o r a g e a r r a y s . The program then goes on t o c a l c u l a t e the p o s i t i o n s and v e l o c i t i e s f o r the next time s t e p . At t h i s p o i n t t h e v e l o c i t i e s i n s t o r a g e a r e f o r time n-1/2 and the c u r r e n t p a r t i c l e v e l o c i t i e s a r e f o r time n+1/2. So t h e n , t o get the v e l o c i t y a t time n, the average i s t a k e n , n n+ 1/ 2 n - 1 / : V = [V + V 3/2 P P P (3.11) T h i s average i s a l s o a c c u r a t e t o second o r d e r . W i t h V p n and the r e s t of the d a t a i n s t o r a g e t h e r e i s a complete s e t of d a t a a t time n t o a n a l y z e . T h i s d a t a however, w i l l be a n a l y z e d d u r i n g the n+1 t i m e s t e p 31 Some Important S c a l i n g Parameters In t h i s s e c t i o n we w i l l d i s c u s s the parameters t h a t p l a y an i m p o r t a n t r o l e i n d e t e r m i n i n g the b e h a v i o u r of the plasma model. They can be c o n s i d e r e d t o be s c a l i n g parameters as they combine c e r t a i n b a s i c plasma c o n s t a n t s and c o n s t a n t s p e r t a i n i n g t o the computer model. We c o n s i d e r now the case of two charge c l o u d s each w i t h a r a d i u s R and a s e p a r a t i o n r . For r£2R t h e Coulomb f o r c e f o r t he c l o u d s i s i d e n t i c a l t o t h a t f o r two p o i n t c h a r g e s . However, f o r r^2R the c l o u d s b e g i n t o o v e r l a p and the e f f e c t i v e p o i n t charge e q u i v a l e n t f o r each c l o u d s t a r t s t o d e c r e a s e . The Coulomb f o r c e i s t h u s m o d i f i e d , d r o p p i n g o f f f o r d e c r e a s i n g s e p a r a t i o n . For r=0 the f o r c e i s z e r o . T h i s i s i n c o n t r a s t t o the case of p o i n t p a r t i c l e s where f o r s m a l l s e p a r a t i o n s the f o r c e becomes l a r g e . I t i s t h e s e s h o r t ( <R ) range i n t e r a c t i o n s where the Coulomb f o r c e i s l a r g e , t h a t account f o r c o l l i s i o n a l e f f e c t s . A consequence then of r e p l a c i n g t h e p o i n t p a r t i c l e w i t h a c l o u d i s t h a t c o l l i s i o n a l e f f e c t s a r e much reduc e d . I t was mentioned i n the d i s c u s s i o n of Debye s c r e e n i n g t h a t the c o l l e c t i v e b e h a v i o u r of a plasma i s determined by l o n g r a n g e (r>X D) Coulomb i n t e r a c t i o n s . Thus i f the c l o u d r a d i u s R 3 ^ the model n e g l e c t s the i n d i v i d u a l a s p e c t s of plasma b e h a v i o u r w h i l e r e t a i n i n g t h e c o l l e c t i v e b e h a v i o u r . I f R i s made much l a r g e r than X p though, the c o l l e c t i v e 32 be h a v i o u r e x h i b i t e d by the model s t a r t s t o d i v e r g e from the p h y s i c a l system ( L a n g d o n ( 7 0 ) ) . For i n t h i s case R r e p l a c e s X D as the boundary between the i n d i v i d u a l and c o l l e c t i v e r egimes. So f o r the CIC scheme, which i s what we used, where R=Ax, the be s t r e s u l t s a r e o b t a i n e d f o r Ax' » X . D As a n a t u r a l consequence of h a v i n g X D = Ax, one of the b a s i c r e q u i r e m e n t s f o r the e x i s t e n c e of a plasma i s f u l f i l l e d , t h a t b e i n g L >> X D, where L i s the plasma l e n g t h . In the model, L = (N - 1 ) A X , and s i n c e X n =* Ax and N u » 1 , the M D M ' plasma l e n g t h w i l l be much l a r g e r than X D. Another s c a l i n g f a c t o r of i n t e r e s t i s d e f i n e d as v M = v t h At/Ax . (3.13) I f vm =1 t h e n a p a r t i c l e w i t h a v e l o c i t y v t h , w i l l t r a v e l a d i s t a n c e of one mesh l e n g t h i n one time s t e p . I f v^ > 1, t r u n c a t i o n e r r o r s a r e i n t r o d u c e d t h r o u g h the d i f f e r e n c i n g e q u a t i o n s . T h i s i s because the v e l o c i t i e s of p a r t i c l e s w i t h v^ £ v M>1 ( v ^ =[v/v t h ] v M ) , a r e u n c o r r e l a t e d t o the e l e c t r i c f i e l d s t he p a r t i c l e s pass t h r o u g h . For i n s t a n c e , i f a p a r t i c l e p o s i t i o n b e f o r e moving i s (j+1/2)Ax and i t s new v e l o c i t y i s v as p a r t l y d e t e r m i n e d by the e l e c t r i c f i e l d a t j and j+1, i t s p o s i t i o n a t the end of the time s t e p w i l l be (j+1/2+v M)Ax. In t h i s t i me s t e p t h e n , t h e p a r t i c l e t r a v e r s e s a c r o s s v c e l l s o b l i v i o u s t o the e l e c t r i c f i e l d s a t the mesh p o i n t s i t pa s s e s by. C l e a r l y then i t i s 33 d e s i r a b l e t o choose the i n i t i a l p arameters so t h a t v M <1. On the o t h e r hand, t h e r e i s no need t o make v M too s m a l l . S i n c e the s p a t i a l r e s o l u t i o n of the f i e l d s i s Ax, i t i s j u s t o v e r k i l l i f v^ i s made any s m a l l e r than 0.1 . The f i n a l s c a l i n g parameter of i n t e r e s t i s w pAt. The f r e q u e n c y of plasma waves i s of the o r d e r of CJ . C o n s e q u e n t l y the s a m p l i n g i n t e r v a l At s h o u l d be < wp ~ 1 i f a c c u r a t e r e s u l t s a r e t o be a c h e i v e d and so we have w pAt < 1. From the p o i n t of view of n u m e r i c a l s t a b i l i t y , i f the system i s t o have a s t a b l e s o l u t i o n , CJ At s h o u l d be l e s s than 2. p For v a l u e s l a r g e r than 2, the s o l u t i o n blows up (see Hockney ( 8 1 ) ) . : To summarize t h e n , t h e p r a c t i c a l e n v e l o p e s of t h e s e t h r e e parameters can be e x p r e s s e d as M P X /Ax - 1 (3.14) 0.1 < v At/Ax < 1 (3.15) co pAt < 1 (3.16) The parameters can a l s o be combined so t h a t ( UpAt )( X D/Ax ) = ( V | h At/Ax ) (3.16) 34 CHAPTER 4 RESULTS LANDAU DAMPING A common t e s t of the a c c u r a c y of a program i n plasma s i m u l a t i o n i s i n how w e l l i t h a n d l e s Landau damping. T h i s s e c t i o n w i l l d e s c r i b e the r e s u l t s of s i m u l a t i o n s t h a t e x h i b i t e d Landau damping. Comparisons of n u m e r i c a l and t h e o r e t i c a l damping v a l u e s w i l l be made. In o r d e r t o s i m u l a t e Landau damping we p e r t u r b the d i s t r i b u t i o n f u n c t i o n f ( r , v , t ) a t time z e r o . So then a t t=0 we have f ( r , v , 0 ) = f o ( r , v , 0 ) + f , ( r , v , 0 ) (4.1) where f o ( r , v , 0 ) i s the u n p e r t u r b e d d i s t r i b u t i o n f u n c t i o n a t t=0 and f , ( r , v , 0 ) i s the p e r t u r b a t i o n . I f the damping i s t o be l i n e a r the p e r t u r b a t i o n must be s m a l l . From t h i s p e r t u r b a t i o n an e l e c t r o s t a t i c o r plasma wave w i l l form. The wave w i l l damp away a t r a t e d e t e r m i n e d by i t s wavelength and the Debye l e n g t h . I t s wavelength i s s e t by the i n i t i a l i z a t i o n p a r a m e t e r s . I f £(x,v , t ) i s the n o r m a l i z e d d i s t r i b u t i o n f u n c t i o n and 35 n ( x , t ) i s the number d e n s i t y of the plasma, f ( x , v , t ) can be w r i t t e n as f ( x , v , t ) = n ( x r t ) ? ( x , v , t ) (4.2) S i n c e the model i s one d i m e n s i o n a l , x has been s u b s t i t u t e d f o r the r v e c t o r and v now r e p r e s e n t s v e l o c i t i e s i n the x d i r e c t i o n o n l y . The f u n c t i o n f ( x , v , t ) can be p e r t u r b e d e i t h e r by p e r t u r b i n g the d e n s i t y n ( x , t ) d i r e c t l y , or by p e r t u r b i n g the v e l o c i t y d i s t r i b u t i o n , f ( x , v , 0 ) . Both t h e s e methods were used i n t h i s work. The d e n s i t y v a r i a t i o n method i s d e s c r i b e d f i r s t . The d e n s i t y n ( x , t ) , i s the e l e c t r o n number d e n s i t y . The i o n number d e n s i t y remains c o n s t a n t throughout the plasma r e g i o n . At the s t a r t of a r u n , t=0, the e l e c t r o n s a r e p l a c e d i n the plasma i n such a way t h a t a s m a l l s i n u s o i d a l v a r i a t i o n i s i n t r o d u c e d i n t o n ( x , 0 ) . The e l e c t r o n d e n s i t y i s g i v e n by n(x,0) = n 0 [ l + a s i n ( k x + <t>) ] 0 < x < L 0 < 4> £ 2TT (4.3) where k i s the wave number of the wave and L i s plasma l e n g t h . F o r l i n e a r damping, a must be s m a l l ( <0.1). As th e c h a rge d e n s i t y i s c a l c u l a t e d o n l y a t mesh p o i n t s , i n t h e program, n(x,0) must be w r i t t e n as 36 n(x,0) = n( ( j - 1 )Ax, 0) = n c [ 1 + a s i n (k ( j - 1 / 2 ) Ax + <f>) ] 1 < j < N M -1 (4.4) where N i s the number of mesh p o i n t s . As w e l l the change has been made from n 0 t o n c . Something t h a t must be kept i n mind when c h o o s i n g the i n i t i a l p a r a m e t e r s i s t h a t because t h e r e i s a f i n i t e number of mesh p o i n t s , t h e number of w avelengths w h i c h t h e mesh can a c c u r a t e l y r e p r e s e n t i s a l s o l i m i t e d . The e l e c r o n s a r e p l a c e d i n the plasma a c c o r d i n g t o E q u a t i o n (4.4) a t the s t a r t of a r u n , t=0. The number of p a r t i c l e s i n t h a t c e l l e q u a l s n ( ( j - 1 ) A x , 0 ) . The s p a c i n g between p a r t i c l e s w i t h i n a c e l l i s u n i f o r m . A l l p a r t i c l e s w i t h i n the plasma a r e g i v e n random v e l o c i t i e s t aken from the M a x w e l l i a n d i s t r i b u t i o n , * f ( x , v , 0 ) = f ( v , 0 ) = [ v t h * e x p ( - v 2 / v t h 2 ) ~ (4.5) The o t h e r method t h a t was used t o g e n e r a t e waves c o n s i s t e d of p e r t u r b i n g ^ ( x , v , 0 ) . To do t h i s , e q u a t i o n (4.5) was a l t e r e d by v a r y i n g s i n u s o i d a l l y i n x, the t h e r m a l v e l o c i t y so, v = v [ 1 + a s i n ( k ( j - l / 2 ) A x + <t>)] th t h 1 <: j < N M-1 a £0.1 (4.6) 37 As mentioned j u s t p r e v i o u s l y , a t the onset of a run the p a r t i c l e s are a s s i g n e d random v e l o c i t i e s t a k e n from the d i s t r i b u t i o n f ( x , v , 0 ) . G i v e n t h e s e c o n d i t i o n s , p a r t i c l e s w i t h i n c e l l j w i l l be g i v e n random v e l o c i t i e s t a k e n from 1 ( ( j - 1 ) A x , v , 0 ) . T h i s w i l l be done f o r a l l N M -1 c e l l s . The p a r t i c l e s a r e p l a c e d i n the plasma w i t h a c o n s t a n t s e p a r a t i o n between p a r t i c l e s . Thus a t t=0 the e l e c t r o n d e n s i t y i s c o n s t a n t t hroughout and n ; = n e . So the charge d e n s i t y i s z e r o everywhere and t h e r e w i l l be no a p p r e c i a b l e waves p r e s e n t . T h i s s i t u a t i o n w i l l not l a s t l o n g . A wave w i t h the same k as i n e q u a t i o n (4.6) w i l l q u i c k l y appear. In e i t h e r c a s e , the wave t h a t i s g e n e r a t e d behaves l i k e a damped o s c i l l a t i o n , d e c a y i n g a t a r a t e d e t e r m i n e d by the the Landau damping c o e f f i c i e n t . From the o utput g e n e r a t e d by v a r i o u s d i a g n o s t i c s , t h e damping r a t e f o r the g e n e r a t e d wave can be c a l c u l a t e d and compared w i t h the t h e o r e t i c a l v a l u e . The method employed i n t h i s work c o n s i s t e d of d e t e r m i n i n g the decay r a t e of the wave i n k space. At r e g u l a r time i n t e r v a l s , t h e F o u r i e r t r a n s f o r m of the e l e c t r i c f i e l d i s t a k e n g i v i n g E ( k ) . A time h i s t o r y of the k spectrum i s t h u s o b t a i n e d . G e n e r a l l y though, a time h i s t o r y of the complete k spectrum i s not r e a l l y r e q u i r e d . I n s t e a d , a few s p e c i f i c k modes, i n p a r t i c u l a r , the p e r t u r b a t i o n mode, a r e chosen a t the s t a r t of a run and f o l l o w e d i n t i m e . The o b v i o u s advantage of t h i s method i s 38 t h a t the e v o l u t i o n of s p e c i f i c k mode can be f o l l o w e d . Damping r a t e s c a l c u l a t e d from averaged q u a n t i t i e s such as the RMS E f i e l d s u f f e r i n t h e i r a c c u r a c y due t o c o n t r i b u t i o n s of s p u r i o u s s i g n a l s from modes o t h e r than the p r i n c i p l e one. T h i s i s not so much a problem f o r ti m e s e a r l y i n a run where the s i g n a l t o n o i s e r a t i o i s s t i l l q u i t e h i g h . For t i m e s t ^ 7 " 1 though, the s i g n a l can become l o s t i n the n o i s e . However, t h e r e i s a n o t h e r s o u r c e of n o i s e t h a t w i l l c o n t r i b u t e t o the k h i s t o r i e s . T h i s source i s a phenomenon known as a l i a s i n g . The minimum wavelength the mesh can r e c o g n i z e i s Ax. The a c c u r a c y of the p a r t i c l e placement though a l l o w s f o r waves w i t h w a velengths l e s s than Ax t o e x i s t i n the plasma. The mesh has no way of r e c o g n i z i n g t h e s e h i g h e r f r e q u e n c y modes. The r e s u l t i s t h a t t h e r e w i l l be waves w i t h k > 27r/Ax t h a t w i l l masquerade as o t h e r modes a d d i n g f a l s e c o n t r i b u t i o n s t o t h e v a r i o u s k h i s t o r i e s . Two t y p i c a l h i s t o r i e s a r e shown i n F i g u r e s 2 and 3. Each shows t h e e v o l u t i o n of t h e p e r t u r b a t i o n mode i n k space. The d e n s i t y a l o n e was i n i t i a l l y p e r t u r b e d i n F i g u r e 2. In F i g u r e 3 the t h e r m a l v e l o c i t y was p e r t u r b e d . As the F o u r i e r t r a n s f o r m r o u t i n e g e n e r a t e s a c u r v e l i k e c o s 2 ( w t ) e " 7 t , t h e p e r i o d of the wave i s the time t a k e n f o r two o s c i l l a t i o n s . I n f i g u r e s 2 and 3 and i n many of the f o l l o w i n g f i g u r e s a number l a b e l l e d as S c a l e i s d i s p l a y e d . I f t h i s number i s p r e s e n t i t means t h a t the v e r t i c a l a x i s 39 Time H i s t o r y Of The P e r t u r b a t i o n In K Space D e n s i t y P e r t u r b e d k/k D = 0.4 N p = 29484 S c a l e = 2.4x10 s 1.0 F i g u r e 2 40 Time H i s t o r y Of The P e r t u r b a t i o n I n K Space Thermal V e l o c i t y P e r t u r b e d k/k = 0.4 N = 29484 D p S c a l e • 7.4x10* F i g u r e 3 41 L e a s t Squares F i t To Peaks In F i g u r e 2 F i g u r e 4 42 has been s c a l e d t o 1.0, so t h a t 1.0 on the v e r t i c a l a x i s i s a c t u a l l y e q u a l t o S c a l e . A f e a t u r e i n F i g u r e 2 t h a t was found t o occur i n a l l s i m u l a t i o n s u s i n g the d e n s i t y p e r t u r b a t i o n scheme i s the anomalously h i g h damping i n g o i n g from the f i r s t peak t o the second peak. T h i s damping i s much h i g h e r than the damping observ e d t o t a k e p l a c e a f t e r the second peak. T h i s i s more c l e a r l y e v i d e n t i n F i g u r e 4 which shows a l e a s t - s q u a r e s f i t of the peak v a l u e s of the o c i l l a t i o n s i n F i g u r e 2. The f i r s t peak s t a n d s w e l l c l e a r of t h e f i t t e d l i n e t h r o u g h t h e o t h e r p o i n t s . T h i s can be a t t r i b u t e d t o the f a c t t h a t i n i t i a l l y , when the p e r t u r b a t i o n i s l o a d e d i n t o the plasma, the p a r t i c l e v e l o c i t i e s a r e u n c o r r e l a t e d w i t h the e l e c t r i c f i e l d s g e n e r a t e d by the p a r t i c l e p l a c e m e n t s . D u r i n g the f i r s t o s c i l l a t i o n the plasma wave and the p a r t i c l e v e l o c i t i e s a d j u s t themselves so t h a t they become s e l f - c o n s i s t e n t . T h i s r e - a d j u s t i n g a c c o u n t s f o r the h i g h e r damping d u r i n g the f i r s t o s c i l l a t i o n . To c a l c u l a t e the damping r a t e s , graphs such as the one F i g u r e 4 were used. The s l o p e of the l e a s t squares f i t l i n e A, d e t e r m i n e s the s c a l e d damping c o e f f i c i e n t 7 /o>p . The r e l a t i o n between the two q u a n t i t i e s i s 7 /wp = A/( 2 - r l o g ( e ) ) . The f i t t i n g was done w i t h o u t u s i n g t h e f i r s t peak v a l u e f o r the reason j u s t mentioned. W h i l e t h e r e i s no e q u i v a l e n t t o e x p e r i m e n t a l e r r o r i n a 43 computer s i m u l a t i o n ( e x c e p t perhaps f o r round o f f e r r o r ) , t h e r e a r e s t a t i s t i c a l f l u c t u a t i o n s t o worry about. T h i s f a c t i s r e f l e c t e d i n t h e f i g u r e s of T a b l e I I I . The runs i n T a b l e I I I were g i v e n an i n i t i a l d e n s i t y p e r t u r b a t i o n so t h a t f = n(a«sin(kx + <f>) ^  ( v ) . The o n l y q u a n t i t i e s a l t e r e d between runs were e i t h e r <j> or the seed. The seed r e f e r s t o the i n i t i a l i z i n g number s u p p l i e d t o the random number g e n e r a t o r , FRAND, a t t h e s t a r t of a r u n . G i v e n the same seed, FRAND w i l l a l w a y s g e n e r a t e the same sequence of numbers. Changing the seed w i l l r e s u l t i n a d i f f e r e n t s e t of random numbers. However, p r o v i d e d t h e r e a r e enough p a r t i c l e s , a M a x w e l l i a n v e l o c i t y d i s t r i b u t i o n s h o u l d be g e n e r a t e d no m a t t e r what the seed i s . As w e l l , s i n c e the plasma i s p e r i o d i c , c h a n g i n g $ s h o u l d not make any r e a l d i f f e r e n c e . W h i l e i n t h e o r y , c h a n g i n g e i t h e r of t h e s e q u a n t i t i e s s h o u l d not r e s u l t i n a change of the r e s u l t s , T a b l e I I I shows t h a t i n p r a c t i c e t h i s i s not the c a s e . C a l c u l a t e d v a l u e s f o r 7/CJ p a r e shown a l o n g w i t h t h e t h e o r e t i c a l v a l u e s . The v a r i a t i o n s w i t h d i f f e r i n g v a l u e s of 0 and t h e seed a r e c l e a r l y e v i d e n t . These f l u c u a t i o n s a r e s t a t i s t i c a l i n n a t u r e due t o the degree of randomness t h a t i s a p a r t of the model. I t i s no s u r p r i s e t h e n , t h a t the w i d e s t d e v i a t i o n s o c c u r f o r the runs where N =10080. When N i s i n c r e a s e d t o p p 29484, t h e r e b y i m p r o v i n g the s t a t i s t i c s , t he d e v i a t i o n s d e c r e a s e and the computed damping r a t e s come c l o s e r t o the T a b l e I I I Landau Damping V a l u e s O b t a i n e d From S i m u l a t i o n s D =0.4 T h e o r e t i c a l v a l u e i s 7/co = 6.613 x 1 7/w px 10" 2 N P SEED 1 0084 29484 50400 5 0.0 4.108 6.534 7.436 10 0.0 8.656 6.864 8.234 18 0.0 7.089 8.015 7.410 30 0.0 8.553 6.420 6.486 0 . 5ir 6.717 6.468 6.835 0 It 4.471 6.494 6.816 0 1 .5rr 3.826 5.845 6.835 0 1 . 8TT 7.760 6.563 7.699 45 the t h e o r e t i c a l v a l u e . ' What i s s u r p r i s i n g i s t h a t when N p i s i n c r e a s e d t o 50000 the numbers do not improve and the d e v i a t i o n s do not d e c r e a s e f u r t h e r . The upshot of a l l t h i s i s t h a t i f t h e damping r a t e i s t o be c a l c u l a t e d f o r a g i v e n w a v e l e n g t h , one run w i l l not s u f f i c e . S e v e r a l runs must be done w i t h d i f f e r e n t seeds or tf>'s and t h e average damping r a t e t a k e n . I n c r e a s i n g t h e number of runs w i l l of c o u r s e improve the f i t , but weighed a g a i n s t t h i s i s the i n c r e a s e i n computer time r e q u i r e d . The modeler must choose where the b a l a n c e l i e s between i n c r e a s i n g a c c u r a c y and i n c r e a s i n g c o s t . In t h i s work, t o determine the damping f o r a p a r t i c u l a r k/k Q, t h e r e would be f o u r runs where the seed=0 and 0 = 0 , 7 r / 2 , IT, 3ir/2. Averages f o r y/up and <Vwp a r e then c a l c u l a t e d from the f o u r r u n s . F i g u r e 5 shows the r e s u l t i n g v a l u e s f o r 7/wp o b t a i n e d , compared a g a i n s t t h e t h e o r e t i c a l c u r v e . F i g u r e 6 does the same f o r u>/cjp . The number of p a r t i c l e s used was the same f o r a l l runs and was e q u a l t o 31500. In both F i g u r e s 5 and 6 t h e d a t a agree q u i t e c l o s e l y w i t h the t h e o r y . In F i g u r e 5 the f i t worsens f o r l a r g e r k/k D (>.4). T h i s i s u n d e r s t a n d a b l e s i n c e as k/k i n c r e a s e s so does the damping. T h i s i n t u r n r e s u l t s i n fewer o s c i l l a t i o n s whose peaks can be used f o r d a t a p o i n t s i n d e t e r m i n i n g y/u , and the a c c u r a c y d e c r e a s e s . 46 F i g u r e 5 47 F i g u r e 6 48 Energy C o n s i d e r a t i o n s For the e l e c t r o s t a t i c system the t o t a l system energy i n a volume V i s g i v e n by The f i r s t term on the r i g h t i s the e l e c t r o s t a t i c energy of the plasma wave, and the second term i s the t o t a l k i n e t i c energy of the e l e c t r o n s i n the system. There i s no i o n c o n t r i b u t i o n t o the k i n e t i c energy s i n c e the i o n s a r e assumed t o be m o t i o n l e s s . I f t h e r e a r e no e x t e r n a l f i e l d s a p p l i e d , the t o t a l system energy s h o u l d remain c o n s t a n t w i t h t i m e . T h i s f a c t was used as a n o t h e r means of t e s t i n g the program. The t o t a l k i n e t i c energy and the e l e c t r o s t a t i c energy were c a l c u l a t e d a t r e g u l a r i n t e r v a l s d u r i n g a r u n . The t o t a l system energy c o u l d then be found by a d d i n g the two q u a n t i t i e s . F i g u r e s 7,8,9 show the time e v o l u t i o n of E K , E £ and E T f o r the same two runs r e f e r r e d t o i n F i g u r e s 2 and 3. The p l o t s of E E and E k c l e a r l y i l l u s t r a t e the o s c i l l a t i o n of E T (4.7) v 49 Energy H i s t o r i e s Thermal V e l o c i t y P e r t u r b e d k/k = 0.4 D 1 F i g u r e 50 Energy H i s t o r i e s D e n s i t y P e r t u r b e d k/k =0.4 D 1 9 7 J, K i n e t i c S c a l e = 1 . 1 x 1 0 7 I i i i 0 1 — i — i — i — i — | — i — i — i — i — i — i — i — i i r 'A, E l e c t r o s t a t i c S c a l e = 3 . 2 x ! 0 5 F i g u r e 8 T o t a l Energy Thermal V e l o c i t y P e r t u r b e d S c a l e = 1 . 1 x 1 0 7 .99985 9997 .9996 F i g u r e 9 52 energy between the d i s t r i b u t i o n f u n c t i o n as r e p r e s e n t e d by the k i n e t i c energy and t h e e l e c t r o s t a t i c energy of the plasma wave. However, d u r i n g b o t h runs t h e t o t a l energy remained r e l a t i v e l y c o n s t a n t . In both runs the a m p l i t u d e of the f l u c u a t i o n s i n t o t a l energy were of the o r d e r of .01% of the average t o t a l energy. T h i s p r o v e d t o be t h e case f o r a l l r u n s . V a l u e s f o r t h e average t o t a l energy, and the s t a n d a r d d e v i a t i o n of E ( t ) / ( E maximum) f o r s e v e r a l runs a r e g i v e n i n T a b l e IV. These were the same runs t h a t were used t o d e t e r m i n e the damping r a t e s i n F i g u r e 5. Each v a l u e i s the average of f o u r r u n s . T a b l e IV Sample E n e r g i e s and D e v i a t i o n s k/k D 0 ( x 1 0 7 ) (E/EU,) average 0.32 1 .206 0.013 0.36 1 .550 0.099 0.40 1 .769 0.008 0.45 1 .570 0.000 0.50 1.411 0.008 0.55 1 .282 0.008 0.60 1 . 174 0.008 53 By k e e p i n g t a b s on the t o t a l energy one has a good way of knowing i f the program i s f u n c t i o n i n g c o r r e c t l y . As an i l l u s t r a t i o n , s e v e r a l runs were made where the o n l y q u a n t i t y v a r i e d was the s i z e of the time s t e p . In c h a p t e r 3 i t was n o t e d the program i s u n s t a b l e f o r w pAt>2. The program w h i l e s t a b l e remains i n a c c u r a t e f o r 1<a) pAt<2. T h i s can be c l e a r l y seen i n F i g u r e 10. Here the t o t a l energy has p l o t t e d i n time f o r s e v e r a l d i f f e r e n t time s t e p s . For w pAt=0.25, the energy remains c o n s t a n t . The time s t e p s i z e used i n most of the Landau damping s i m u l a t i o n s i n t h i s work used time s t e p s of a p p r o x i m a t e l y t h i s s i z e . As the time s t e p s i z e i n c r e a s e s the t o t a l energy s t a r t s t o show more v a r i a t i o n and an o v e r a l l i n c r e a s e i n v a l u e . F i n a l l y , f o r a>pAt > 2 the system blows up as e v i d e n c e d by the r a p i d and l a r g e i n c r e a s e i n t o t a l energy. R e v e r s i b i l i t y So l o n g as c o l l i s i o n a l e f f e c t s a r e n e g l i g i b l e Landau damping i s a r e v e r s i b l e p r o c e s s . Thus i f a t some time t 0 the system i s r e v e r s e d ( i . e . f o r a l l p a r t i c l e s v -> -v a t time t 0 ) the plasma s h o u l d r e t u r n t o i t s o r i g i n a l s t a t e a t time 2 t 0 . As t h e program models a c o l l i s i o n l e s s plasma i t s h o u l d a l s o e x h i b i t r e v e r s i b i l i t y . F i g u r e 11 shows the r e s u l t s of F i g u r e 1 0 55 one such s i m u l a t i o n . I n t h i s c a s e , the d e n s i t y was p e r t u r b e d w i t h the parameters chosen so t h a t k/k D=0.4 . The s i m u l a t i o n proceeded i n a normal f a s h i o n ' u n t i l t=5r p ( T p = the plasma p e r i o d ) . At t h i s time a l l v e l o c i t i e s were r e v e r s e d . The s i m u l a t i o n was then c o n t i n u e d t o time t = 1 6 r p . F i g u r e 11 i s a h i s t o r y of the p e r t u r b a t i o n i n k space. From t=0 t o t = l O r the p l o t c l e a r l y m i r r o r s about t=5r p the r e v e r s a l p o i n t . The damping r a t e from t=0 t o 5 r p was c a l c u l a t e d and was found t o be e q u a l t o y/up=-0.1482. The growth r a t e f o r t=5r p t o t = l 0 r p was d e t e r m i n e d t o be 7/a>p=+0.1498. The peak a t t = l 0 r p came t o w i t h i n 97% of the i n i t i a l peak v a l u e a t t=0. A f t e r the plasma has r e t u r n e d t o i t s i n i t i a l s t a t e , the wave once a g a i n s t a r t s t o damp away. An i n t e r e s t i n g f e a t u r e of F i g u r e 11 i s t h a t b o t h the peak a t t=0 and the peak a t t = l 0 r p s t a n d w e l l c l e a r of the n e i g h b o u r i n g peaks. The h i g h damping r a t e e x h i b i t e d . d u r i n g the f i r s t o s c i l l a t i o n i s m i r r o r e d by a h i g h growth r a t e a t t = l 0 r . p A s i m i l a r run was made where the t h e r m a l v e l o c i t y was p e r t u r b e d . T h i s run a l s o d i s p l a y e d r e v e r s i b i l i t y though not as w e l l as the p r e v i o u s c a s e . Time H i s t o r y Of The P e r t u r b a t i o n I n K Space D e n s i t y P e r t u r b e d k/k =0.4 N a 31500 V e l o c i t i e s Reversed A t t= 5T P in CTN 57 TWO STREAM INSTABILITY A n o t h e r t e s t of the program t h a t was c a r r i e d out was t o see how w e l l i t s i m u l a t e d the two stream i n s t a b i l i t y . To do t h i s , the . e l e c t r o n s were i n i t i a l i z e d so t h a t the plasma c o n s i s t e d of two c o u n t e r s t r e a m i n g e l e c t r o n beams of e q u a l d e n s i t y ( n 0 / 2 ) and e q u a l and o p p o s i t e v e l o c i t i e s v 0 i n a f i x e d i o n background. The s u s c e p t i b i l i t y of t h i s phenomenon t o l o n g w a v e l e n g t h i n s t a b i l i t y i s c l e a r l y shown i n F i g u r e 1 2 ( a - d ) . F i g u r e s 12 and 13 a r e phase space p l o t s of random sample of p a r t i c l e s . V e l o c i t y i s p l o t t e d i n v t h u n i t s and p o s i t i o n i s p l o t t e d i n terms of t h e plasma l e n g t h L. Each dot r e p r e s e n t s one p a r t i c l e . I n t h i s s i m u l a t i o n t h e r e was no p e r t u r b a t i o n so t h a t f o ( v ) = \ «(v - v o ) + ^_ 6(v + v ) ( 4 . 8 ) 2 . 2 where v 0=9.7v t h . W h i l e t h e r e was no a c t u a l wave p r e s e n t a t t=0, k can be g i v e n by k=2ir/L where L i s t h e plasma l e n g t h . For t h e parameters chosen i n F i g u r e 12, kv o=0.97. That k=2ir/L i s c l e a r from F i g u r e s I 3 ( a - b ) . Here v 0 = 11v j h so t h a t f o r t h e same parameters used i n the p r e v i o u s run k v 0 =1.1 which i s above the i n s t a b i l i t y t h r e s h o l d (=1.0). T h i s e x p l a i n s why F i g u r e s 13(a-b) remain r e l a t i v e l y unchanged. F i g u r e s 14 and 15 show p l o t s f o r t h e t o t a l k i n e t i c and 58 e l e c t r i c , and t o t a l system energy f o r kv o=0.97 s i m u l a t i o n . The l a r g e f i r s t o s c i l l a t i o n i n F i g u r e 14 marks the onset of the i n s t a b i l i t y . As i n the p r e v i o u s s e c t i o n , the t o t a l system energy s h o u l d remain c o n s t a n t w i t h t i m e . W h i l e t h e r e was an i n c r e a s e i n t o t a l energy f o r kv o=0.97 the i n c r e a s e was s m a l l ( < 2% ). P l o t s of t h e e n e r g i e s f o r kv 0=1.1 were not i n c l u d e d as t h e r e was o n l y s m a l l v a r i a t i o n s f o r E , E , or E K . A s e r i e s of s i m u l a t i o n s were done where k v 0 v a r i e d from 0.1 t o 0.9. For the s e runs X=L, a=0.1 and k=0.1 . S t a r t i n g a t 1.0, v 0 was in c r e m e n t e d i n s t e p s of ol.0v t h. I n i t i a l v e l o c i t y d i s t r i b u t i o n s were once a g a i n monoenergetic (v=±v 0). H i s t o r i e s of the k i n e t i c and e l e c t r i c energy f o r each run were o b t a i n e d . I n p l o t t e d form t h e y a r e s i m i l a r i n form t o F i g u r e 14. They do d i f f e r from each o t h e r however i n the time t h a t the l a r g e f i r s t o s c i l l a t i o n a p p e a r s . T a b l e V shows the v a l u e s o b t a i n e d f o r 7 ' f o r the v a r i o u s v a l u e s of k v 0 , where 7 ' i s the i n v e r s e of the time ( i n r p " 1 u n i t s ) t h a t the peak of the l a r g e f i r s t o s c i l l a t i o n o c c u r s . 20, Phase Space P l o t For Two Stream I n s t a b i l i t y k v Q = 0.9 7 t = QTp F i g u r e 12(a-b) 60 20 o A •20 * 1 , . . - - , , x . „ . t= 4 0. .5 1. 20 t = 6 r p o H Ma** * l • • .<;r- .. . • • • H i > . • 20, 0. .5 F i g u r e 12(c-d) 61 Phase Space P l o t For Two Stream I n s t a b i l i t y k v 0 = 1.1 t = O . O r P 20.1 > oH t = 6 . 0 T P 201 > o H •20. 0 . —r .5 X F i g u r e 13(a-b) K i n e t i c And E l e c t r o s t a t i c Energy For Two Stream I n s t a b i l i t y kv 0 = 0.97 S c a l e = 8.1X10 6 F i g u r e 14 63 T o t a l Energy For Two Stream I n s t a b i l i t y kv = 0.97 S c a l e = 8.2X10 6 0 4 8 t / T p F i g u r e 15 64 T a b l e V Computed Growth Rates For Two Stream I n s t a b i l i t y k v 0 7' 0.1 0.4129 0.2 0.4546 0.3 0.5359 0.4 0.5552 0.5 0.5598 0.6 0.5549 0.7 0.5476 0.8 0.5073 0.9 0.4343 Comparing T a b l e V and T a b l e I I we see t h a t 7 ' and 7 ( t h e o r y ) both peak i n the neighbourhood of k v 0 ^ 0.6. In the next s i m u l a t i o n t h e program was t e s t e d f o r r e v e r s i b i l t y . The plasma was p e r t u r b e d so t h a t v Q = V q ( 1 + a s i n ( 2 7 r x / x ) ) (4.9) where X=L, a=0.1 and v 0=6v so t h a t kv o=0.6 . F u r t h e r m o r e , f o r t h i s case the v e l o c i t y d i s t r i b u t i o n c o n s i s t e d of two m a x w e l l i a n s c e n t e r e d about ±v 0 i n s t e a d of two monoenergetic 65 beams a t ±v 0. At t = 5 r p a l l v e l o c i t i e s were r e v e r s e d and the s i m u l a t i o n c o n t i n u e d t o t = l O r p . F i g u r e s 16 and 17 show the r e s u l t i n g h i s t o r i e s a l l t h r e e e n e r g i e s . I f the energy v a l u e s a t the v e r y b e g i n n i n g and end of the run a r e d i s r e g a r d e d , the p l o t s m i r r o r v e r y w e l l about t=5T p, the r e v e r s a l p o i n t . Good r e v e r s i b i l i t y i s shown. The h i g h e r growth r a t e f o r kv o=0.6 compared t o the growth r a t e f o r kv o=0.97. i s e v i d e n c e d i n F i g u r e 17 by the e a r l i e r appearance of the i n s t a b i l i t y . In F i g u r e 14 where kv o=0.97, the f i r s t peak o c c u r s a t t ^ 2.5T p w h i l e i n F i g u r e 17 the f i r s t peak o c c u r s a t t * 1.7T P . 66 K i n e t i c And E l e c t r o s t a t i c Energy For Two Stream I n s t a b i l i t y V e l o c i t i e s R e v e r s e d a t t-5r p kv 0 = 0.60 S c a l e = 5.1X10 7 K i n e t i c 0.4 J F i g u r e 16 67 T o t a l Energy For Two Stream I n s t a b i l i t y V e l o c i t i e s Reversed a t t=5r p kv B 0.60 S c a l e = 5.2x10' 0 1.0 J 0.9 8 ]—i—i 1 1 1 1 1 j — i 1 1 1 1 1 1 1—i 1 1—j 0 4 8 F i g u r e 17 68 CHAPTER 5 CONCLUSIONS In t h i s c o n c l u d i n g c h a p t e r a r e v i e w of what was a c c o m p l i s h e d i s g i v e n . As w e l l , a few t h o u g h t s on f u t u r e a p p l i c a t i o n s of t h i s type of program are o f f e r e d . The g o a l of t h i s work has been t o d e v e l o p a w o r k i n g e l e c t r o s t a t i c CIC code and by u s i n g v a r i o u s means, det e r m i n e i t s p h y s i c a l r e l i a b i l i t y . The program was used t o s i m u l a t e two w e l l known plasma phenomena; Landau damping and the two stream i n s t a b i l i t y . In the case of Landau damping, damping r a t e s f o r s e v e r a l v a l u e s of k/k D were c a l c u l a t e d from program o u t p u t . These were i n good agreement w i t h the t h e o r e t i c a l v a l u e s . The f r e q u e n c i e s w of t h e plasma waves d e t e r m i n e d by t h e program were a l s o found t o be c l o s e t o t h e o r y . For the two stream i n s t a b i l i t y the program agreed w i t h the t h e o r e t i c a l t h r e s h o l d q u i t e w e l l , showing no growth f o r k v 0 > 1 . S i n c e t h e s e a r e b o t h r e v e r s i b l e p r o c e s s e s , the program was a l s o t e s t e d f o r r e v e r s i b i l i t y . The s i m u l a t i o n s e x h i b i t e d t h i s f e a t u r e q u i t e c l e a r l y f o r b o t h phenomena. F i n a l l y , as d E ( t o t a l ) / d t = 0 i n e i t h e r phenomenon, checks were made f o r the c o n s t a n c y of t h e t o t a l energy. In Landau 69 damping E was seen t o be c o n s t a n t w i t h i n .01%. For the two stream i n s t a b i l i t y the t o t a l energy v a r i a t i o n was l a r g e r ( <2%). In e i t h e r case the d e v i a t i o n s were i n the form of f l u c t u a t i o n s and d i d not r e p r e s e n t any r e a l growth i n the t o t a l energy. T h i s p r o v i d e d good e v i d e n c e of t h e s t a b l e n a t u r e of the program. The program t h e n , or one l i k e i t , w i l l be a u s e f u l t o o l i n f u t u r e plasma s t u d i e s . W h i l e i n i t s p r e s e n t s t a t e i t i s l i m i t e d t o s i m u l a t i n g e l e c t r o s t a t i c phenomenon w i t h f i x e d i o n s , i t c o u l d e a s i l y be adapted t o l o o k a t a wide range of plasma p r o c e s s e s . By r e p l a c i n g t h e c o n s t a n t p o s i t i v e background charge w i t h moving p o s i t i v e p a r t i c l e s , s i m u l a t i o n s c o u l d be done of phenomenon f o r w h i c h t h e i o n m o t i o n i s i m p o r t a n t . T h i s would r e l a t i v e l y s i m p l e t o implement. However, t h i s e n t a i l s c a l c u l a t i n g p o s i t i o n s and v e l o c i t i e s f o r 2N p a r t i c l e s and computing time i s b a s i c a l l y d o u b l e d . The time s c a l e i s of o r d e r u and runs w i t h t h e s e time d u r a t i o n s can become r a t h e r e x p e n s i v e as t h e time s t e p s i z e i s s t i l l l i m i t e d by the c o n d i t i o n <o pAt< 1 . By a l t e r i n g t h e method of s o l v i n g f o r E t h e s t e p s i z e c o n s t r a i n t c o u l d be r e l a x e d . One such method i s t o s o l v e f o r E v i a an i m p l i c i t s o l v i n g r o u t i n e (Mason ( 8 1 ) ) . Use of t h i s method a l l o w s f o r w p A t » 1 . Another a r e a of study t h a t would be of i n t e r e s t i s l a s e r - p l a s m a i n t e r a c t i o n s . To do t h e s e type of s i m u l a t i o n s e l e c t r o - m a g n e t i c f i e l d s must be f i r s t i n c o r p o r a t e d i n t o the 70 code. As t h i s would e n t a i l s o l v i n g M a x well's e q u a t i o n s i n d i f f e r e n c e form though, a new l i m i t known as the Courant c o n d i t i o n i s i n t r o d u c e d l i m i t i n g the s i z e of the t ime s t e p . T h i s c o n d i t i o n e x i s t s due t o t h e n u m e r i c a l b e h a v i o u r of f i n i t e d i f f e r e n c e s o l u t i o n s of p a r t i a l d i f f e r e n t i a l e q u a t i o n s . In o r d e r f o r t h e EM PIC code t o be s t a b l e , the s t e p s i z e i s l i m i t e d t o At < Ax/2C. I f the a p p l i e d f i e l d s a r e h i g h enough, r e l a t i v i s t i c c o r r e c t i o n s t o the v e l o c i t y c a l c u l a t i o n s would a l s o have t o be made. As a f i n a l remark i t s h o u l d be n o t e d t h a t s i n c e the number of p a r t i c l e s g r e a t l y exceeds the number of mesh p o i n t s , most of the computing time used i n a run w i l l be spent on moving the p a r t i c l e s and f i n d i n g t h e i r v e l o c i t i e s . S p e c i a l c a r e t h e n , must be taken t o o p t i m i z e the p a r t i c l e moving r o u t i n e . I t was f o r t h i s r e a s o n t h a t t h e p a r t i c l e moving r o u t i n e i n t h i s work was w r i t t e n i n Assembler. On the average run t h i s reduced computing time by about 10% when compared t o a s i m i l a r run u s i n g an a l l F o r t r a n code. 71 REFERENCES B i r d s a l l , C.K., and F u s s , D., J o u r n a l of C o m p u t a t i o n a l P h y s i c s 3, 494 ( 1 9 6 9 ) . B u r g e r , P., J o u r n a l of A p p l i e d P h y s i c s 3 £ , 1938 ( 1 9 6 5 ) . Canosa, J . , J o u r n a l of C o m p u t a t i o n a l P h y s i c s J_3, 158 ( 1 9 7 3 ) . Chen, F.,F., I n t r o d u c t i o n To Plasma P h y s i c s , Plenum P r e s s , New York, (1974) Dawson, J.M., P h y s i c s of F l u i d s 5, 158 ( 1 9 6 2 ) . Dawson, J.M., Methods i n C o m p u t a t i o n a l P h y s i c s 9 , 1 ( 1 9 7 0 ) . Dawson, J.M., Reviews of Modern P h y s i c s 5 5 , 403 ( 1 9 8 3 ) . Hockney, R.W., P h y s i c s of F l u i d s 9 , 1826 ( 1 9 6 6 ) . Hockney, R.W., and Eastwood, J.W., Computer S i m u l a t i o n s U s i n g P a r t i c l e s , M c G r a w - H i l l , New York, ( 1 9 8 1 ) . I c h i m a r u , S., B a s i c P r i n c i p l e s Of Plasma P h y s i c s ( F r o n t i e r s i n P h y s i c s ) , W.A. Benjamin I n c . , Reading M a s s e c h u s e t t s , ( 1 9 7 3 ) . Landau, L., J o u r n a l of P h y s i c s (USSR) J_0, 25 ( 1 9 4 6 ) . Langdon, A.B. and L a s i n s k i , B.F., Methods i n C o m p u t a t i o n a l P h y s i c s J_6, 327 ( 1976) . Langdon, A.B. and B i r d s a l l , C.K., P h y s i c s of F l u i d s J_3, 21 15 ( 1 9 7 0 ) . Mason, R.J., J o u r n a l of C o m p u t a t i o n a l P h y s i c s 4 1 , 233 ( 1 9 8 1 ) . 72 APPENDIX A DETERMINING THE TOTAL ENERGY When c a l c u l a t i n g the system energy the d i f f e r e n c e between the p h y s i c a l s i t u a t i o n and model environment must be taken i n t o a c c o u n t . In t h i s s e c t i o n a d e s c r i p t i o n of how the e n e r g i e s were d e t e r m i n e d i s g i v e n We r e c a l l from e q u a t i o n (4.7) t h a t the t o t a l energy i n a volume V f o r the e l e c t r o s t a t i c system i s g i v e n by r 1.2 1 . E T 7 • 8ir 2 ' : 1 v where N i s the number of e l e c t r o n s i n the volume V. In the case of a one d i m e n s i o n a l plasma model however, because the plasma i s i n f i n i t e i n the y-z p l a n e i t makes more sense t o c a l c u l a t e the energy s u r f a c e d e n s i t y . W i t h t h i s i n mind t h e e l e c t r o s t a t i c energy becomes I E 2 (x) I ff 1 dxdydz = E _ l / d y d z (A.2) 8ir where E f f i s t h e e l e c t r o s t a t i c energy s u r f a c e d e n s i t y [ e r g s / c m 2 ] . To c a l c u l a t e t h e ^ k i n e t i c energy we f i r s t note t h e number of p a r t i c l e s i n a s i m u l a t i o n N p , i s not e q u a l t o N. 73 T h i s i s because each s i m u l a t i o n p a r t i c l e r e p r e s e n t s an i n f i n i t e ( i n y-z) sheet of charge w i t h a s u r f a c e number d e n s i t y ff . S i n c e the model (which c o n t a i n s N p p a r t i c l e s over a d i s t a n c e L) r e p r e s e n t s a plasma w i t h number d e n s i t y n 0 , we have N a P n = n 0 (A.3) N Q can be d e t e r m i n e d from the i n p u t parameters and so a v a l u e f o r <rn can be o b t a i n e d . T h e r e f o r e a model p a r t i c l e w i t h a v e l o c i t y v c o n t r i b u t e s mea^v2/2 t o the k i n e t i c energy s u r f a c e d e n s i t y . So i f E t o . i s t h e t o t a l energy s u r f a c e d e n s i t y , we have N p = Ea + ! > ° n ^ 2 v i 2 [ e r g s / c m 2 ] 2 i = 1 (A.4) T h i s i s what i s c a l c u l a t e d i n the program. APPENDIX B NAMELIST PARAMETERS A l l t he parameters f o l l o w i n g a r e i n the n a m e l i s t I n i t i a l i z i n g Parameters INPSWOO) = Chooses type of i n i t i a l c o n d i t i o n s . = 1 D e n s i t y p e r t u r b e d = 2 Thermal v e l o c i t y p e r t u r b e d = 3 Two stream case A1 = a m p l i t u d e of i n i t i a l p e r u r b a t i o n A2 = a m p l i t u d e of second p e r t u r b a t i o n ANE1 = number of wavelengths i n plasma f o r i n i t i a l p e r t u r b a t i o n ANE2 = number of wavelengths i n plasma f o r second p e r t u r b a t i o n PH11 = phase a n g l e ( i n p i u n i t s ) f o r i n i t i a l p e r t u r b a t i o n PHI2 = phase a n g l e ( i n p i u n i t s ) f o r second p e r t u r b a t i o n VO = Streaming v e l o c i t y f o r the two stream case MONOEN = Used i n c o n j u c t i o n w i t h t h e two stream c a s e . I f MONOEN=1 the v e l o c i t y d i s t r i b u t i o n w i l l c o n s i s t of two streams a t v=+/- VO. I f MONOEN=0 the v e l o c i t y d i s t r i b u t i o n w i l l c o n s i s t of two m a x w e l l i a n s c e n t e r e d about +/- VO SEED = P r o v i d e s the i n i t i a l i z i n g seed f o r FRAND. NP = Number of p a r t i c l e s UTHIN = Thermal v e l o c i t y ( i n u n i t s of c) OMEGIN = Plasma frequency ( i n u n i t s of 3.55x10«•15 ANL = Le n g t h of plasma ( i n u n i t s of 5.3x10*«-5 (cm)) MESH = Number of mesh p o i n t s i n g r i d CFACT = Determines s i z e of t h e time s t e p , DT DT=DX/(CFACT«Vth) ICON = T e l l s the program i f t h i s run i s a c o n t i n u a t i o n of a p r e v i o u s run 1=yes, i t ' s a c o n t i n u a t i o n 0=no NT = Number of time s t e p s f o r the run 75 NTPULS = I f NTPULS i s g r e a t e r than 0, when the program reaches the time step=NTPULS a second p e r t u r b a t i o n i s i n t r o d u c e d . I f NTPULS <-2, when the program r e a c h e s the time step=-NTPULS a l l v e l o c i t i e s a r e r e v e r s e d . O t h e r w i s e NTPULS s h o u l d =-2. INPSW(12) = S e t s the type of p a r t i c l e boundary c o n d i t i o n s used. 1 p a r t i c l e bounces back o f f edge w i t h a new random v e l o c i t y 2 p a r t i c l e bounces back o f f edge w i t h v e l o c i t y r e v e r s e d 3 p a r t i c l e i s r e i n t r o d u c e d a t o t h e r s i d e w i t h a new random v e l o c i t y i n the d i r e c t i o n i t was moving 4 p a r t i c l e i s r e i n t r o d u c e d a t o t h e r s i d e w i t h v e l o c i t y unchanged. Program uses Assembler v e r s i o n of p a r t i c l e mover. Program uses F o r t r a n v e r s i o n of p a r t i c l e mover. Program uses b o t h on a l t e r n a t e time s t e p s . MOVRTN = 1 2 3 D i a g n o s t i c Parameters F o u r i e r D i a g n o s t i c s INVLFT = F o u r i e r T r a n s f o r m of E ( x ) ta k e n e v e r y INVLFT time s t e p NFTSEL = Number of F o u r i e r t r a n s f o r m h i s t o r i e s t h a t a r e t o be r e c o r d e d . I f NFTSEL=-1 the complete spectrum w i l l be r e c o r d e d . MFTSEL(I) = C o n t a i n s which K's a r e t o be r e c o r d e d . I f NFTSEL=-1 the program f i l l s i n MFTSEL by i t s e l f . F i e l d H i s t o r y INPSW(11) = 0 o n l y h i s t o r i e s of c e r t a i n mesh p o i n t s t a k e n 1 o n l y rms v a l u e s f o r E and RHO w i l l be c a l c u l a t e d . 2 b o t h of the above w i l l be done. INVLFL = F i e l d h i s o r y t a k e n e v e r y INVLFT time s t e p NMSEL = Number of mesh p o i n t s t o be r e c o r d e d 76 MSHSEL(I) = C o n t a i n s which mesh p o i n t s a r e t o be r e c o r d e d . I f NFTSEL=-1 the program f i l l s i n MFTSEL by i t s e l f . Energy D i a g n o s t i c s INVLEN = Energy ( k i n e t i c , e l e c t r o s t a t i c and t o t a l ) c a l c u l a t e d e v e r y INVLEN time s t e p . Phase Space P l o t s INVLXV = Phase space p l o t g e n e r a t e d f o r e v e r y INVLXV time s t e p . NSAMP = T h i s i s the number of p a r t i c l e s used i n the phase space p l o t VMAX2 = Range of v e l o c i t y a x i s on p l o t i s +/- vmax2. Other P l o t s INVLPL = P l o t s done f o r e v e r y INVLPL time s t e p . INPSW(1) = 2 P l o t f ( v ) 1 asks u s e r f i r s t 0 Don't p l o t f ( v ) VMAX = f ( v ) p l o t t e d from 0 t o VMAX ( i n V t h u n i t s ) INPSW(2) = 2 P l o t E ( x ) 1 ask u s e r f i r s t 0 Don't p l o t E ( x ) T r a c e r s INVLTR = V e l o c i t y and p o s i t i o n of t r a c e r p a r t i c l e s r e c o r d e d e v e r y INVLTR time s t e p NUMTR = Number of t r a c e r s TRP(I) = T r a c e r p o s i t i o n s ( i n u n i t s of t h e plasma l e n g t h L) TRV(I) = T r a c e r v e l o c i t i e s ( i n V t h u n i t s ) )s C C C PROGRAM NAME:PICSTAT C C WRITTEN BY ARNE ALFHEIM, 1983 C C THIS IS A PARTICLE IN CELL (PIC) CODE. C FOR A REVIEW OF PIC CODES SEE DAWSON, J.,M., C REVIEWS OF MODERN PHYSICS [ 5 5 ] , 403 1983. C C THE PROGRAM SIMULATES A ONE DIMENSIONAL PLASMA 77 C WITH FIXED IONS. IT IS ALSO AN ELECTRO-STATIC C SIMULATION AS THERE ARE NO APPLIED FIELDS. C c======================================================== c C THE FOLLOWING IS A DESCRIPTION OF SOME OF THE C COMMONLY USED VARIABLES. THESE VARIABLES DO NOT C OCCUR IN THE NAMELIST. C c c C ARRAYS C C UX CONTAINS THE PARTICLE VELOCITIES C XP CONTAINS THE PARTICLE POSITIONS C C EX CONTAINS THE ELECTRIC FIELD AT THE MESH POINTS C RHO CONTAINS THE CHARGE DENSITY AT THE MESH POINTS C C RHOD IS A DUMMY ARRAY USED FOR TEMPORARY STORAGE C C THE NEXT THREE ARRAYS ARE USED IN THE DIAGNOSTICS. C C VOLD IS USED TO STORE THE VELOCITIES FOR THE C PREVIOUS TIME STEP C XOLD IS USED TO STORE THE POSITIONS FOR THE C PREVIOUS TIME STEP C RHOLD IS USED TO STORE THE CHARGE DENSITIES FOR THE C PREVIOUS TIME STEP C c c C CONSTANTS (BASED IN CGS UNITS) C C ALAMDA= LASER WAVELENGTH ( HALF MICRON) C THIS SERVES AS A BASE UNIT FOR THE PLASMA C LENGTH. C OMEGA = LASER FREQUENCY. SERVES AS BASE UNIT C FOR THE PLASMA FREQUENCY C TIME = REAL TIME (SECONDS) C C XL = PLASMA LENGTH C DX = CELL LENGTH C DT = TIME STEP LENGTH C MESH = NUMBER OF MESH POINTS C MESHM1= MESH-1 (NUMBER OF CELLS) C C OMEGAP= PLASMA FREQUENCY (1/SEC) C UTH = THERMAL VELOCITY (CM/SEC) C CID = NUMBER OF PARTICLE PER CELL C DENSUR= n«XL/Np where n i s the plasma d e n s i t y C C EOVERM= ELECTRON CHARGE TO MASS RATIO 78 C QET EOVERM'DT C QDX e - n / c i d C where e i s the e l e c t r o n charge C C = VELOCITY OF LIGHT C PI = 3.14 C C C=========================================== c COMMON /PARAR/ XP(100000), UX(100000) COMMON /FLDAR/ EX(1000), RHO(1000), RHOD(lOOO) COMMON /MSHCON/ MESH, MESHM1, NP, DX, XL, DT COMMON /EMCON/ EOVERM, QET, QDX, C, PI COMMON /LASCON/ OMEGA, ALAMDA, TIME COMMON /DENCON/ OMEGAP, UTH, CID, DENSUR COMMON /DGNS/ VMAX, INPSW(30), MSHSEL(200), NMSEL, 1 MFTSEL(129), NFTSEL, INVLPL, INVLFT, INVLFL, 2 INVLEN, INVLXV COMMON /PHASPL/ VMAX2, VMINPH, VMAXPH, I SAMP, NSAMP COMMON /PARC/ A1, ANE1, PHI 1, A2, ANE2, PHI 2, 1 NTPULS, NUMTR, TRP(30), TRV(30), INVLTR, 2 UTHIN, OMEGIN, ANL, VO COMMON /STRGE/ VOLD(100000),XOLD(100000),RHOLD(1000) NAMELIST /INDAT/ ICON, NT, UTHIN, OMEGIN, CFACT,ANL, 1 NP, MESH, A1 , ANE1, PHI 1 , A2, ANE2, PHI 2, 2 NTPULS, NUMTR, TRP, TRV, INPSW, MFTSEL, 3 NFTSEL, MSHSEL, NMSEL, INVLPL, INVLTR, 4 INVLFL, INVLFT, INVLEN, VMAX, SEED, VO, 5 VMAX2, INVLXV, NSAMP, MOVRTN, MONOEN C C C THE FOLLOWING ARE DEFAULT VALUES FOR THE C INPUT PARAMETERS FROM THE NAMELIST INDAT C MOVRTN = 1 MONOEN = 0 ICON = 0 SEED = 0.0 OMEGIN = .2 CFACT = 4.0 UTHIN = .01 ANL = 2.0 NP = 6000 MESH = 128 A1 = . 1 PHI 1 = .5 ANE1 = 10. A2 = .1 ANE2 = 1. PHI 2 = .5 NTPULS = -2 NT = 10 NUMTR = 0 NMSEL = 1 NFTSEL = 0 VMAX = 10. VMAX2 = 0. VO = 0. NSAMP = 1000 INVLTR = 10000 INVLPL = 10000 INVLFL = 10000 INVLFT = 10000 INVLXV = 10000 INPSW(1) = 0 INPSW(2) = 0 INPSW(3) = 0 INPSWO0) = 1 INPSW(11) = 1 INPSW(12) = 4 READ IN INITIALIZING DATA NTSC = 0 READ (5,1NDAT) INITIALIZE RANDOM NUMBER GENERATOR XXX = RAND(SEED) IF (NTPULS .GE. 0) NT = NT - 1 NTS = -1 IF (ICON .EQ. 1) NTS = 0 IF THIS IS A CONTINUATION THE E FIELDS FOR THE FIRST STEP OF THE RUN ALREADY EXIST, SO NTS=0 . IF THIS IS NOT A CONTINUATION THE INITIAL FIELDS HAVE TO BE DETERMINED WITHOUT INTERFERING WITH THE INITIAL PARTICLE DISTRIBUTION. THIS IS DONE BY SETTING NTS=-1 FOR THIS FIRST STEP THE PARTICLES ARE NOT MOVED NOR IS TIME ADVANCED. ONLY RHO AND E ARE CALCULATED. SUBROUTINE INVAL CALCULATES VARIOUS CONSTANTS THAT WILL BE USED IN THE RUN SUCH AS DT, DX, AND CID CALL INVAL(TIMEC, ICON, NTSC, CFACT) IF THIS IS A NOT A CONTINUATION (ICON=0), TIMEC=0. OTHERWISE TIMEC IS THE TIME AT WHICH THE PREVIOUS RUN STOPPED. 80 TIME = TIMEC IF (ICON .EQ. 0) GO TO 10 CALL DATA(NTS) C C DATA CONTAINS VALUES FOR EX,XP,UX FOR THE LAST C TIME STEP OF THE PREVIOUS RUN. C GO TO 20 10 CONTINUE C C MYZERO ZEROS ARRAYS OF DIMENSION, MESH C CALL MYZERO(EX, MESH) C C THE NEXT STEP IS TO SET UP THE INITIAL DISTRIBUTION C DEPENDING ON THE VALUE OF INPSW(10), ONE OF THE C INIT ROUTINES IS CALLED. INITIAL PARTICLE C VELOCITIES AND POSITIONS ARE ASSIGNED HERE. C I F THE RUN IS A CONTINUATION, THIS SECTION IS C BYPASSED. C IF (INPSWO0) .EQ. 1) CALL INIT1 IF (INPSWO0) .EQ. 2) CALL INIT2 IF (INPSWdO) .EQ. 3) CALL INIT3 (MONOEN) CALL TRACE C 20 CONTINUE C c C SOME DIAGNOSTIC QUANTITIES ARE SET UP IN INDGNS. C CALL INDGNS(ICON) C C C C C THE LINES THAT FOLLOW, COMPRISE THE MAIN LOOP, C STEPPING FORWARD IN TIME THE PARTICLES. C 30 CALL MYZERO(RHO, MESH) 40 CONTINUE C C IF (INPSW(10).NE.3) GO TO 45 IF ((NTS .EQ. IABS(NTPULS)).AND.(NTPULS .LT. - 2)) 1 CALL REVERS(NP) C 50 CONTINUE C C SUBROUTINE PULSE GENERATES A SECOND DENSITY C PERTURBATION AT A TIME STEP NTPULS. C FOR STEPS WHERE NT.NE.NTPULS, PULSE DOES NOTHING C CALL PULSE(NTS) 81 C C C THE HEART OF THE CODE IS IN STEPON. C STEPON CONTAINS THE PARTICLE MOVER C AND THE ELECTRIC FIELD SOLVER. C CALL STEPON(NTS, NTPULS, NTSC, MOVRTN) C C THIS IS JUST TO SAY IT'S WORKING C C C C C C C C C C C C C C-C C C C C IF (M0D(NTS,1) .EQ. 0) WRITE (6,60) NTS IF (NTS .LT. NT) GO TO 30 C C C c C THE PROGRAM HAS RUN ITS COURSE AND NOW WRITES INTO C A FILE THE DATA NECESSARY TO RESUME RUNNING AT SOME C OTHER (USER'S) TIME. AS WELL ANY OTHER DIAGNOSTICS C DESIRED CAN BE DONE AT THIS POINT. C C C C C I F THERE WAS A SECOND PULSE (NTPULS>0) THEN NTS C MUST BE CORRECTED. IF (NTPULS .GE. 0) NTS = NTS - 1 NTS = NTS + NTSC C OUTPUT WRITES INTO DEVICE 1 THE DATA REQUIRED FOR C CONTINUATION. CALL OUTPUT(NTS) STOP 60 FORMAT (16) END SUBROUTINE STEPON(NTS, NTPULS, NTSC, MOVRTN) COMMON /PARAR/ XP(100000), UX(100000) COMMON /FLDAR/ EX(lOOO), RHO(1000), RHOD(lOOO) COMMON /LASCON/ OMEGA, ALAMDA, TIME COMMON /MSHCON/ MESH, MESHM1, NP, DX, XL, DT BEGIN = 1. IF ((NTS.EQ.- 1) .OR. (NTS.EQ.NTPULS)) BEGIN=0. IF (BEGIN .EQ. 1.) GO TO 10 IF THIS IS AN INITIAL TIME STEP OR IF THERE IS A SECOND PULSE IN THIS STEP THEN THE ELECTRIC FIELD MUST BE FIRST SET UP BEFORE MOVING THE PARTICLES. THIS IS WHAT NOMOVE DOES. CALL NOMOVE GO TO 30 CONTINUE MOVEPA IS THE ASSEMBLER VERSION OF THE PARTICLE MOVER MOVEPF IS THE FORTRAN VERSION OF THE PARTICLE MOVER THE PARTICLE MOVER CALCULATES THE NEW PARTICLE POSITIONS AND VELOCITIES AND THE CHARGE DENSITY. IF (MOVRTN .EQ. 1) CALL MOVEPA IF (MOVRTN .EQ. 2) CALL MOVEPF IF (MOVRTN .LT. 3) GO TO 30 IF MOVRTN EQUAL 3 PROGRAM ALTERNATELY CALLS MOVEPF AND MOVEPA THIS IS USEFUL IN COMPARING THE RELATIVE SPEEDS OF MOVEPF AND MOVEPA USING TIMETALLY IN DEBUG IF (MOD(NTS,2) .EQ. 0) GO TO 20 CALL MOVEPF GO TO 30 CONTINUE CALL MOVEPA CONTINUE SWITCH COPIES RHO INTO RHOD CALL SWITCH(RHO, RHOD) TIME = TIME + DT • BEGIN ARRCEN CENTERS RHO TO THE MIDDLE OF THE CELL CALL ARRCEN(RHO) IF (NTS .LT. 0) GO TO 40 83 C VELCEN AVERAGES VOLD TO TIME NTS C CALL VELCEN(NTS, NP) C C PREDGN DOES ALL THE DIAGNOSTICS C CALL PREDGN(NTS, NTSC) C C ARROLD SAVES EX(N) (IN EOLD) AND V(N+1/2) C (INTO VOLD) C 40 CALL ARROLD(NTS, NP) C C ESTAT SOLVES FOR EX(N+1) USING POISSON EQUATION C CALL ESTAT NTS = NTS + 1 RETURN END C C C c C MOVEPF CALCULATES XP(N+1),UX(N+1),RHO(N+1), C SUBROUTINE MOVEPF C COMMON /PARAR/ XP(100000), UX(100000) COMMON /FLDAR/ EX(1000), RHOO000), RHOD( 1000) COMMON /MSHCON/ MESH, MESHM1, NP, DX, XL, DT COMMON /EMCON/ EOVERM, QET, QDX, C, PI COMMON /DENCON/ OMEGAP, UTH, CID, DENSUR C CALCULATE INTERPOLATION CONSTANTS FOR X AT LEVEL N. C RHE =0.0 C --DO FOR ALL THE PARTICLES DO 60 IE = 1, NP C --LOCATE THE CELL C XTEST = X P ( I E ) / DX J = INT(XTEST) + 1 PX = XTEST + 1.0 - FLOAT(J) PX1 = 1.0 - PX C C THE INTERPOLATED FIELDS AT THE PARTICLE. C J1 = J + 1 EXI = PX1 • EX( J ) + PX • EX(J1) C EXI=EX(j)/2.+EX(J1)/2. UX(IE) = UX(IE) + QET • EXI XP( I E ) = X P ( I E ) + UX(IE) • DT C C BOUNDARY CONDITIONS. 84 C C IF PARTICLE GOES PAST EITHER EDGE ROUTINE EDGE C IS CALLED. C C RIGHT-HAND BOUNDARY. C 10 IF (XP(IE) .LE. XL) GO TO 20 CALL EDGE(XP(IE), U X ( I E ) , XL, UTH, -1.0) C C C LEFT-HAND BOUNDARY. C 20 IF (XP(IE) .GE. 0.0) GO TO 30 CALL EDGE(XP(IE), U X ( I E ) , XL, UTH, 1.0) C C 30 CONTINUE C C CALCULATE NEW INTERPOLATION CONTANTS. C 40 XTEST = X P ( I E ) / DX J = INT(XTEST) + 1 PX = XTEST + 1.0 - FLOAT(J) PX1 = 1.0 - PX J1 = J + 1 C c c RHO(J) = RHO(J) + PX1 • QDX IF ( J .EQ. MESH) GO TO 50 RHO (J1 ) = RHO (J1 ) + PX • QDX 50 CONTINUE C C FOR PERIODIC BOUNDARY CONDITIONS RHO(1)=RHO(MESH) C C THE NEXT LINES ENSURE THIS WILL BE SO. IF ( ( J .EQ. 1) .OR. ( J .EQ. MESH)) RHE = RHE + 1 PX1 • QDX IF ( J .EQ. MESHM1) RHE = RHE + PX • QDX C 60 CONTINUE C RHO(1) = RHE RHO(MESH) = RHE C C CION IS THE ION CHARGE DENSITY C CION = QDX • CID C DO 70 J = 1, MESH RHO(J) = RHO(J) - CION 70 CONTINUE C RETURN END C C c c c c SUBROUTINE EDGE(X, V, XL, UTH, FSIGN) C THIS ROUTINE HANDLES PARTICLES THAT HAVE C GONE PAST THE PLASMA EDGES. C INPSW(12)=1 BOUNCE BACK BOUNDARY CONDITIONS C VELOCITY IS REFLECTED WITH NEW C RANDOM VELOCITY C INPSW(12)=2 BOUNCE BACK BOUNDARY CONDITIONS C VELOCITY IS REFLECTED C INPSW(12)=3 PERIODIC BOUNDARY CONDITIONS C VELOCITY IS GIVEN NEW RANDOM C VALUE BUT DIRECTION IS UNCHANGED C INPSW(12)=4 PERIODIC BOUNDARY CONDITIONS C VELOCITY UNCHANGED C c REAL•4 X, V, XL, UTH, FSIGN COMMON /DGNS/ VMAX, INPSW(30), MSHSEL(200), NMSEL, 1 MFTSEL(129), NFTSEL, INVLPL, INVLFT,INVLFL, 2 INVLEN, INVLXV COMMON /PARC/ A1, ANE1, PHI 1, A2, ANE2, PHI 2, 1 NTPULS, NUMTR, TRP(30), TRV(30), INVLTR, 2 UTHIN, OMEGIN, ANL, VO XL2 = 2. • XL XADD = 0. IF (FSIGN .EQ. - 1.) XADD = XL2 IF (INPSW(12) .GT. 2) GO TO 10 C C THE FOLLOWING LINES SIMULATE BOUNCE BACK C BOUNDARY CONDITIONS C X = XADD - X IF (INPSW(12) .EQ. 1) GO TO 20 C V = -V GO TO 30 C C C THE FOLLOWING LINES SIMULATE PERIODIC C BOUNDARY CONDITIONS C 10 X = FSIGN • XL + X IF (INPSW0 2) .EQ. 4) GO TO 30 86 20 V = -FSIGN • ABS(RDMV(UTH)) IF (INPSW(10) .EQ. 3) V = V - FSIGN • VO • UTH 30 RETURN END C C c c C AVERAGES ARRAY VALUES TO MIDPOINTS C C SUBROUTINE ARRCEN(ARRAV) DIMENSION ARRAV(1) COMMON /MSHCON/ MESH, MESHM1, NP, DX, XL, DT C DO 10 I = 1, MESHM1 ARRAV(I) = (ARRAV(I) + ARRAV(I + 1)) • 0.5 10 CONTINUE C RETURN END C C C C C SOLVES FOR THE NEW ELECTRIC FIELD VIA POISSONS C EQUATION. C SUBROUTINE ESTAT COMMON /FLDAR/ EX(lOOO), RHOO000), RHOD(1000) COMMON /MSHCON/ MESH, MESHM1, NP, DX, XL, DT EX(1) = 0.0 ECORSM = 0.0 C DO 10 J = 1, MESHM1 E X ( J + 1) = E X ( J ) + RHO(J) • DX ECORSM = ECORSM + E X ( J + 1) 10 CONTINUE C C ECOR IS A CORRECTION TERM. IDEALLY IT SHOULD C BE ZERO. ECOR = ECORSM / FLOAT(MESH) C DO 20 J = 1, MESH EX ( J ) = E X ( J ) - ECOR 20 CONTINUE C RETURN END C C C 87 C C TRACE IS CALLED IF THERE ARE TO BE TRACER PARTICLES C TRACE LOOKS FOR PARTICLES THAT ARE CLOSEST TO THE C POSITIONS REQUESTED IN TRP. THEN THE PARTICLE C IS ASSIGNED THE VELOCITY IN TRV. THESE PARTICLES C ARE THEN SWITCHED WITH THE PARTICLES AT THE END OF C UX AND XP ARRAYS. SUBROUTINE TRACE DIMENSION ITRSOO) COMMON /PARAR/ XP(100000), UX(100000) C COMMON /MSHCON/ MESH, MESHM1, NP, DX, XL, DT COMMON /DENCON/ OMEGAP, UTH, CID, DENSUR COMMON /PARC/ A1, ANE1, PHI 1, A2, ANE2, PHI 2, 1 NTPULS, NUMTR, TRP(30), TRV(30), INVLTR, 2 UTHIN, OMEGIN, ANL, VO C NPP = NP - NUMTR + 1 IF (NUMTR .EQ. 0) RETURN C DO 10 J = 1, NUMTR ITRS(J ) = 0 10 CONTINUE C c C TRP AND TRV MUST FIRST BE SWITCHED TO CGS UNITS C DO 20 J = 1, NUMTR TRP(J) = TRP(J) • XL TRV(J) = TRV(J) • UTH 20 CONTINUE . C DO 100 J = 1, NUMTR IR = 0 C 30 CONTINUE C C SCAN FOR THE PARTICLE IT IS NEAREST TO. DO 40 I = 1, NP IF (TRP(J) .LT. X P ( D ) GO TO 50 40 CONTINUE C C THIS LINE IS IN CASE POSITIONS ARE ENTERED C INCORRECTLY TRP(J) = XL • .5 GO TO 30 50 AG = TRP(J) - X P ( I - 1) AL = X P ( I ) - TRP(J) C C WHICH PARTICLE IS CLOSER7 C X P ( I ) OR X P ( I - I ) ( X P ( I ) > T R P ( J ) ) I F (AG .LT. AL) GO TO 60 IT = I 88 C C C C IMP = -1 GO TO 70 60 IT = I - 1 IMP = 1 70 DO 80 JT = 1, NUMTR IF THIS PARTICLE HAS ALREADY BEEN PICKED FOR SWITCHING, TAKE THE NEXT CLOSEST ONE. IF (ITRS(JT) .EQ. IT) GO TO 90 80 CONTINUE C c c c c c-c c c c c c c c ITRS(J) = IT GO TO 100 90 IR = IR + 1 IMP = -1 • IMP IT = IT - IMP -GO TO 70 100 CONTINUE IR SWITCH THE TRACER PARTICLES WITH THE ONES AT THE END OF THE ARRAYS. DO 110 J = 1, NUMTR IPT = NP + J - NUMTR XT = XP(IPT) ITR = ITRS(J) XP(IPT) = XP(ITR) XP(ITR) = XT UX(IPT) = TRV(J) 110 CONTINUE 120 RETURN END AT THE END OF A RUN SUBROUTINE DATA WRITES ONTO DEVICE 2 THE DATA REQUIRED TO CONTINUE THE RUN IF SO DESIRED. IF THE RUN IS A CONTINUATION OF SOME PREVIOUS RUN, DATA READS FROM DEVICE 1 THE DATA REQUIRED FOR RESTARTING. SUBROUTINE DATA(NTS) COMMON /PARAR/ XP(100000), UX(100000) COMMON /FLDAR/ EX( 1000), RHOOOOO), RHOD( 1000) COMMON /MSHCON/ MESH, MESHM1, NP, DX, XL, DT COMMON /EMCON/ EOVERM, QET, QDX, C, PI COMMON /LASCON/ OMEGA, ALAMDA, TIME COMMON /DENCON/ OMEGAP, UTH, CID, DENSUR COMMON /DGNS/ VMAX, INPSW(30), MSHSEL(200), NMSEL, 89 C C 1 MFTSEL(129), NFTSEL, INVLPL, INVLFT,INVLFL, 2 INVLEN, INVLXV READ (1) (EX(I),1=1,MESH) 10 READ (1) (XP(J),J=1,NP) READ (1) (UX(J),J=1,NP) 20 GO TO 50 ENTRY OUTPUT(NTS) OMEGAP = OMEGAP / OMEGA UTH = UTH / C ANLAM = XL / ALAMDA WRITE (2,70) DT, OMEGAP, UTH, ANLAM, NP, MESH WRITE (2,60) TIME, NTS WRITE (2) (EX(I),1=1,MESH) 30 WRITE (2) (XP(J),J=1,NP) WRITE (2) (UX(J),J=1,NP) 40 CONTINUE 50 RETURN 60 FORMAT (G16.8, 110) 70 FORMAT (4G16.8, 2110) END C C PULSE GENERATES A SECOND PULSE LIKE IN INIT1 AT C TIMESTEP NTPULS PROVIDED THAT NTPULS>1 C c SUBROUTINE PULSE(NTS) C REAL* 4 LAMP, PHI, CHMOOOO), X( 1000) INTEGER-4 XAXT /' X V , YAXT /' XS '/ COMMON /PARAR/ XP(100000), UX(100000) COMMON /MSHCON/ MESH, MESHM1, NP, DX, XL, DT COMMON /EMCON/ EOVERM, QET, QDX, C, PI COMMON /LASCON/ OMEGA, ALAMDA, TIME COMMON /DENCON/ OMEGAP, UTH, CID, DENSUR COMMON /PARC/ A1, ANE1, PHI 1, A2, ANE2, PHI 2, 1 NTPULS, NUMTR, TRP(30), TRV(30), INVLTR, 2 UTHIN, OMEGIN, ANL, VO IF (NTS .NE. NTPULS) RETURN IP = 0 PHI 2 = PHI 2 • PI C DO 20 MN = 1, MESHM1 I = MN - 1 SN = SIN((FLOAT(I)«DX + DX/2.)•ANE2•2.-Pl/XL+PHI2) CED = CID • ( 1 . + A2«SN) TLOP = CED - FLOAT(INT(CED)) LOP = 0 IF (TLOP .GE. 0.5) LOP = 1 NPM = INT(CED) + LOP CHM(MN) = FLOAT(NPM) - CID 90 XSPM = DX / FLOAT(NPM) C DO 10 J - 1, NPM IP = IP + 1 IF ( I P .GT. NP) GO TO 30 J J = J - 1 XTEMP = FLOAT(JJ) • XSPM X P ( I P ) = XTEMP + FLOAT(I) • DX 10 CONTINUE C 20 CONTINUE C GO TO 40 30 LOP = 10 40 IDIF = NP - IP IF (IDIF .LE. 0) GO TO 60 C DO 50 I = 1, IDIF I P = IP + 1 XP ( I P ) = FRAND(0.) • XL 50 CONTINUE C 60 CONTINUE C DO 70 K = 1, MESHM1 70 X(K) = FLOAT(K - 1) / FLOAT(MESHM1) C CALL GPLOT(CHM, X, MESHM1, XAXT, YAXT, 5., 0.) 80 CONTINUE 90 RETURN END C C C c C VELCEN OBTAINS VELOCITY AT TIME N C IT ALSO STORES THE OLD VELOCITY THAT WILL BE USED C IN THE NEXT TIME STEP TO CALCULATE V(N) C IT ONLY DOES THIS WHEN THE DIAGNOSTICS DESIRED C REQUIRE IT C ^ SUBROUTINE VELCEN(NTS, NP) COMMON /PARC/ A1, ANE1, PHI 1, A2, ANE2, PHI 2, 1 NTPULS, NUMTR, TRP(30), TRV(30), INVLTR, 2 UTHIN, OMEGIN, ANL, VO COMMON /STRGE/ VOLD(100000),XOLD(100000),RHOLD(1000) COMMON /DGNS/ VMAX, INPSW(30), MSHSEL(200), NMSEL, 1 MFTSEL0 2 9 ) , NFTSEL, INVLPL, INVLFT, INVLFL, 2 INVLEN, INVLXV COMMON /PARAR/ XP(100000), UX(100000) COMMON /FLDAR/ E X ( I O O O ) , RHO(lOOO), RHOD(1000) C C 91 IF (MOD(NTS,INVLEN) .EQ. 0) GO TO 40 IF (MOD(NTS,INVLXV) .EQ. 0) GO TO 40 IF (MOD(NTS,INVLTR) .NE. 0) GO TO 80 NPP = NP - NUMTR + 1 C IF (NTS .NE. IABS(NTPULS)) GO TO 20 IF (NTPULS .EQ. - 2) GO TO 20 C DO 10 I = NPP, NP VOLD(I) = UX(I) 10 CONTINUE C GO TO 80 20 CONTINUE C DO 30 I = NPP, NP VOLD(I) = (VOLD(I) + U X ( I ) ) • 0.5 30 CONTINUE C GO TO 80 C 40 IF (NTS .NE. IABS(NTPULS)) GO TO 60 IF (NTPULS .EQ. - 2) GO TO 60 C DO 50 I = 1, NP VOLD(I) = UX(I) 50 CONTINUE C GO TO 80 60 CONTINUE C DO 70 I = 1, NP VOLD(I) = (VOLD(I) + U X ( D ) / 2. 70 CONTINUE C 80 RETURN C c c C ARROLD SAVES THE DATA REQUIRED FOR DIAGNOSTICS C IN THE NEXT TIME STEP. C ENTRY ARROLD(NTS,NP) C NTP = NTS + 1 IF (MOD(NTP,INVLFT) .EQ. 0) CALL SWITCH(RHOD, RHOLD) IF (MOD(NTP,INVLFL) .EQ. 0) CALL SWITCH(RHOD, RHOLD) IF (MOD(NTP,INVLEN) .EQ. 0) GO TO 100 IF (MOD(NTP,INVLXV) .EQ. 0) GOTO 100 IF (MOD(NTP,INVLTR) .NE. 0) GO TO 130 NPP = NP - NUMTR + 1 C DO 90 I = NPP, NP 92 C C C C c c c VOLD(I) = UX(I) XOLD(I) = X P ( I ) 90 CONTINUE GO TO 130 100 DO 110 I = 1, NP VOLD(I) = UX(I) 110 CONTINUE IF (MOD(NTP,INVLXV) .NE. 0) GO TO 130 DO 120 I = 1, NP XOLD(I) = X P ( I ) 120 CONTINUE 130 RETURN END C C C NOMOVE IS USED IN THE FIRST TIME STEP TO CALCULATE C THE ELECTRIC FIELDS OF THE INITIAL DISTRIBUTION. C THE PARTICLES ARE NOT MOVED SUBROUTINE NOMOVE C COMMON /PARAR/ XP(100000), UX(100000) COMMON /FLDAR/ EX(lOOO), RHOO000), RHOD(lOOO) COMMON /MSHCON/ MESH, MESHM1, NP, DX, XL, DT COMMON /EMCON/ EOVERM, QET, QDX, C, PI COMMON /DENCON/ OMEGAP, UTH, CID, DENSUR C CALCULATE INTERPOLATION CONSTANTS FOR X C RHE = 0.0 C --DO FOR ALL THE PARTICLES DO 20 IE = 1, NP C XTEST = X P ( I E ) / DX J = INT(XTEST) + 1 PX = XTEST + 1.0 - FLOAT(J) PX1 = 1.0 - PX RHO(J) = RHO(J) + PX1 • QDX IF ( J .EQ. MESH) GO TO 50 J1=J+1 RH0(J1) = RH0(J1) + PX • QDX 50 CONTINUE C c c c I F ( ( J .EQ. 1) .OR. ( J .EQ. MESH)) RHE « RHE + 93 1 PX1 • QDX IF ( J .EQ. MESHM1) RHE = RHE + PX • QDX 20 CONTINUE C RHO(1) = RHE RHO(MESH) = RHE CION = QDX • CID C DO 30 J = 1, MESH RHO(J) = RHO(J) - CION 30 CONTINUE C c RETURN END C C C C INVAL CALCULATES SOME OF THE CONSTANTS SUBROUTINE INVAL(TIMEC, ICON, NTSC, CFACT) COMMON /MSHCON/ MESH, MESHM1, NP, DX, XL, DT COMMON /EMCON/ EOVERM, QET, QDX, C, PI COMMON /LASCON/ OMEGA, ALAMDA, TIME COMMON /DENCON/ OMEGAP, UTH, CID, DENSUR COMMON /PARC/ A1, ANE1, PH11, A2, ANE2, PHI 2, 1 NTPULS, NUMTR, TRP(30), TRV(30), INVLTR, 2 UTHIN, OMEGIN, ANL, VO TIMEC = 0. PI = 4. • ATAN(1.0) C = 3.0E+10 ALAMDA = 0.53E-04 OMEGA = 2 . • PI • C / ALAMDA IF (ICON .EQ. 0) GO TO 10 READ (1,80) CFACT, OMEGIN, UTHIN, ANL, NP, MESH WRITE (6,30) WRITE (6,80) CFACT, OMEGIN, UTHIN, ANL, NP, MESH READ (1,40) TIMEC, NTSC 10 CONTINUE UTH = UTHIN • C MESHM1 = MESH - 1 XL = ANL • ALAMDA DX = XL / FLOAT(MESHM1) DENR1D = FLOAT(NP) / XL DENSR = DENR1D •• 3 OMEGAP = OMEGIN • OMEGA DT = DX / (CFACT*UTH) EMASS = 9.10E-28 ECH = 4.803250E-10 EOVERM = ECH / EMASS DENS = (OMEGAP*•2) / (4.*PI•EOVERM*ECH) DENSUR = DENS • XL / FLOAT(NP) CID = FLOAT(NP) / FLOAT(MESHM1) QET = EOVERM • DT QDX = (OMEGAP••2) / (CID^EOVERM) DEBYEP = UTH / OMEGAP WRITE (6,50) DX, DEBYEP C C 20 WRITE (6,60) WRITE (6,70) DX, XL, DENS, DENSUR, CID, QDX RETURN 30 FORMAT (/'DT, OMEGAP, UTH, ANLAM, NP, MESH') 40 FORMAT (G16.8, 110) 50 FORMAT (/'DX = ',G16.8, 5X, 'DEBYE LENGTH ', G16.8) 60 FORMAT (/'DX, XL,DENS, DENSUR, CED QDX') 70 FORMAT (3G16.8, /, 3G16.8) 80 FORMAT (4G16.8, 2110) END C C c c c C INIT1 SETS UP THE INITIAL DISTRIBUTION TO SIMULATE C LANDAU DAMPING. IN THIS CASE THE ELECTRON DENSITY C VARIES SINUSIODALLY WITH X, CENTERED ABOUT CID. C THE ELECTRON DENSITY IS CONSTANT INSIDE A CELL. C THE VELOCITIES ARE TAKEN RANDOMLY FROM A NORMAL C MAXWELLIAN DISTRIBUTION. C SUBROUTINE INIT1 C REAL•4 LAMP, PHI, CHM(1000), X(1000) INTEGER-4 XAXT / ' X */, YAXT /'• XS '/ COMMON /PARAR/ XP(100000), UX(100000) COMMON /MSHCON/ MESH, MESHM1, NP, DX, XL, DT COMMON /EMCON/ EOVERM, QET, QDX, C, PI COMMON /LASCON/ OMEGA, ALAMDA, TIME COMMON /DENCON/ OMEGAP, UTH, CID, DENSUR COMMON /PARC/ A1, ANE1, PHI 1, A2, ANE2, PHI 2, 1 NTPULS, NUMTR, TRP(30), TRV(30), INVLTR, 2 UTHIN, OMEGIN, ANL, VO 10 I P = 0 PHI = PHI 1 • PI PI2 = PI • 2. XLM = XL - DX C DO 30 MN = 1 , MESHM1 I = MN - 1 SN = SIN((FLOAT(I)•DX + DX/2.)'ANE1«PI2/XL+PHI) CED = CID • ( 1 . + A1»SN) C C CHECK TO SEE HOW LARGE A FRACTION OF PARTICLE IS C LEFTOVER. I F .GE.0.5 THEN NPM (THE NUMBER OF C ELECTRONS PER MESH) IS INCREMENTED BY 1 95 TLOP = CED - FLOAT(INT(CED)) LOP = 0 IF (TLOP .GE. 0.5) LOP = 1 NPM = INT(CED) + LOP C C ARRAY CHM RECORDS NUMBER OF ELECTRON IN EXCESS OF C ION NUMBER (CID) PER CELL CHM(MN) = FLOAT(NPM) - CID XSPM = DX / FLOAT(NPM + 1) C DO 20 J = 1 , NPM IP = IP + 1 IF ( I P .GT. NP) GO TO 40 XTEMP = FLOAT(J) • XSPM XP( I P ) = XTEMP + FLOAT(I) • DX 20 CONTINUE C 30 CONTINUE C GO TO 50 40 LOP = 10 50 IDIF = NP - IP IF (IDIF .LE. 0) GO TO 70 C DO 60 I = 1, IDIF IP = IP + 1 XP ( I P ) = FRAND(0.) • XL 60 CONTINUE C 70 WRITE (6,120) I P , LOP, I D I F , MN C C ASSIGN INITIAL VELOCITIES C DO 80 I = 1, NP UX(I) = RDMV(UTH) 80 CONTINUE C DO 90 K = 1, MESHM1 90 X(K) = FLOAT(K - 1) / FLOAT(MESHM1) C C A PLOT OF THE INITIAL PERTURBATION IS PLOTTED CALL GPLOT(CHM, X, MESHM1, XAXT, YAXT, 5., 0.) 100 CONTINUE WRITE (6,110) RETURN C c c c 110 FORMAT ( I X , 'INIT1') 120 FORMAT (4110) END INIT2 SETS UP THE INITIAL DISTRIBUTION TO SIMULATE LANDAU DAMPING. IN THIS CASE THE THERMAL VELOCITY VALUE VTHM, VARIES SINUSIODALLY IN X, CENTERED ABOUT UTH. VTHM IS USED WHEN THE RANDOM VELOCITIES ARE GENERATED. THE ELECTRON DENSITY IS CONSTANT OVER THE WHOLE PLASMA LENGTH. SUBROUTINE INIT2 REAL'4 LAMP, PHI, 'CHM(lOOO), X(1000) COMMON /PARAR/ XP(100000), UX(100000) COMMON /MSHCON/ MESH, MESHM1, NP, DX, XL, DT COMMON /EMCON/ EOVERM, QET, QDX, C, PI COMMON /LASCON/ OMEGA, ALAMDA, TIME COMMON /DENCON/ OMEGAP, UTH, CID, DENSUR COMMON /PARC/ A1, ANE1, PHI 1, A2, ANE2, PHI 2, 1 NTPULS, NUMTR, TRP(30), TRV(30), INVLTR, 2 UTHIN, OMEGIN, ANL, VO 0 NPM = INT(CID) XSP = DX / CID PHI = PHI1 • PI PI 2 = PI • 2. XLM = XL - DX I = 0 DO 30 MN = 1, MESHM1 J = MN - 1 SN = SIN((FLOAT(J)-DX + DX/2.)•ANE1•PI2/XL+PHI) VTHM = UTH • ( 1 . + A1•SN) XTEMP = DX • FLOAT(J) DO 20 I I = 1, NPM 1 = 1 + 1 IF (I ,GT. NP) GO TO 40 X P ( I ) = XTEMP + FLOAT(II - 1) • XSP UX(I) = RDMV(VTHM) 0 CONTINUE 0 CONTINUE GO TO 50 IF SOMETHING GOES WRONG AND THE ROUTINE TRIES TO INITIALIZE MORE PARTICLES THAN NP THEN IT COMES HERE AND WRITES PARTICLE AND MESH NUMBER. 0 WRITE (6,60) I I , MN 0 CONTINUE WRITE (6,70) C RETURN 60 FORMAT ( ' PARTICLE tr EXCEEDED,II AND MESHir = ', 216) 70 FORMAT ( I X , 'INIT2') END C C c • c c C INIT3 SETS UP THE PARTICLES FOR THE TWO STREAM C INSTABILITY. THE DENSITY IS CONSTANT EVERYWHERE. C IF MONOEN=1 THE VELOCITY DISTRIBUTION CONSISTS C OF TWO MONOENERGETIC STREAMS AT +/- VO. C IF MONOEN=0 THE VELOCITY DISTRIBUTION CONSISTS C OF TWO STREAMS AT +/- VO EACH WITH A MAXWELLIAN C SPREAD OF UTH. C IF THERE IS A PERTURBATION (A1>0) VO VARIES C SINUSIODALLY IN X. C C SUBROUTINE INIT3(MONOEN) C REAL• 4 LAMP, PHI, CHM(lOOO), XdOOO) COMMON /PARAR/ XP(100000), UX(100000) COMMON /MSHCON/ MESH, MESHM1, NP, DX, XL, DT COMMON /EMCON/ EOVERM, QET, QDX, C, PI COMMON /LASCON/ OMEGA, ALAMDA, TIME COMMON /DENCON/ OMEGAP, UTH, CID, DENSUR COMMON /PARC/ A1, ANE1, PHI 1, A2, ANE2, PHI 2, 1 NTPULS, NUMTR, TRP(30), TRV(30), INVLTR, 2 UTHIN, OMEGIN, ANL, VO NPM = INT(CID) XSP = DX / CID PHI = PHI1 • PI PI 2 = PI • 2. XLM = XL - DX UTHS = UTH IF (MONOEN .EQ. 1) UTHS =0.0 1 = 0 C DO 20 MN = 1, MESHM1 J = MN - 1 SN = SIN((FLOAT(J)«DX + DX/2.)•ANE1•PI2/XL+PHI) VOP = VO • ( 1 . + A1«SN) XTEMP = DX • FLOAT(J) C DO 10 I I - 1, NPM 1 = 1 + 1 X P ( I ) = XTEMP + FLOAT(II - 1) • XSP 98 C SINCE THE SAMPLE PARTICLES USED IN THE PHASE SPACE C PLOT ARE PICKED AT A REGULAR INTERVAL, TO GET A C GOOD STATISTICAL SAMPLE THE PARTICLE IS ASSIGNED C RANDOMLY +/- VO SW = FRAND( 0 . CO SIGN = 1. IF (SW .LT. 0.5) SIGN = - 1 . U X ( l ) = RDMV(UTHS) + VOP • UTH • SIGN 10 CONTINUE C 20 CONTINUE C C WRITE (6,60) C RETURN 60 FORMAT (1X, 'INIT3') END C C C C C C SUBROUTINE REVERS(NP) COMMON /PARAR/ XP(100000), UX(100000) C DO 10 I = 1, NP UX(I) = -UX(I) 10 CONTINUE C WRITE (6,20) RETURN 20 FORMAT (1X, 'REVERSE') END C SUBROUTINE MYZERO(A, NDIM) DIMENSION A O ) DO 10 J = 1, NDIM A ( J ) = 0 .0 10 CONTINUE RETURN END C SUBROUTINE SWITCH(ARRIN, ARROUT) REAL•4 ARRIN(1), ARROUT(1) COMMON /MSHCON/ MESH, MESHM1, NP, DX, XL, DT 99 C DO 10 I = 1, MESH ARROUT(I) = ARRIN(I) 10 CONTINUE C RETURN END C C C C C C FUNCTION RDMV GENERATES RANDOM VELOCITIES C FROM A GAUSSIAN DISTRIBUTION WITH ZERO MEAN C AND A STANDARD DEVIATION EQUAL TO THE THERMAL C VELOCITY. C C FUNCTION RDMV(UTH) Z = FRAND(0.) C SIGN=1. C IF (Z .LT.0.5) SIGN=-1. C Z=ABS(2.•(Z -.5)) C RDMV=SIGN-UTH-SQRT(-ALOG(Z)) C THETA=FRAND(0.)'6.283 C RDMV=RDMV•COS(THETA) RDMV = UTH • AGAU(Z) RETURN END C 100 c c C DIAGNOSTICS PACKAGE C c C INDGNS SETS UP THE NECESSARY ARRAYS C AND VARIABLES C C PREDGN DOES ALL THE ACTUAL DIAGNOSTICS C C C SUBROUTINE INDGNS(ICON) IMPLICIT REAL-4(A - H,0 - Z,X) COMMON /DENCON/ OMEGAP, UTH, CID, DENSUR COMMON /DGNS/ VMAX, INPSW(30), MSHSEL(200), NMSEL, 1 MFTSEL0 2 9 ) , NFTSEL, INVLPL, INVLFT, 2 INVLFL, INVLEN, INVLXV COMMON /PHASPL/ VMAX2,VMINPH,VMAXPH,I SAMP,NSAMP COMMON /PARC/ A1, ANE1, PHI 1, A2, ANE2, PHI 2, 1 NTPULS, NUMTR, TRP(30), TRV(30), INVLTR, 2 UTHIN, OMEGIN, ANL, VO COMMON /LASCON/ OMEGA, ALAMBA, TIME COMMON /PLCOM/ X1(1000), X2(1000), X 3 O 0 0 0 ) , AXL1, AXL2, AXL3 COMMON /MSHCON/ MESH, MESHM1, NP, DX, XL, DT C C C C C PLOTS SECTION, INITIALIZATION C C XXHH,XKKH ARE HORIZONTAL AXES ARRAYS C INPSW IS USED AS A SWITCHING ARRAY FOR VARIOUS C OUTPUT DECISIONS, GRAPHS ETC. C C C X1 X ARRAY FOR FIELD PLOTS IN X SPACE C X2 X ARRAY FOR FOURIER TRANSFORMS OF THE FIELDS. C X3 X ARRAY FOR VELOCITY DISTRIBUTION DO 10 I = 1, MESH X 1 ( I ) = FLOAT(I - 1) / FLOAT(MESHM1) 10 CONTINUE C M22 = MESH / 2 C DO 20 I = 1, M22 X 2 ( I ) = FLOAT(I - 1) 20 CONTINUE C DO 30 I = 1, 51 X 3 ( I ) = (VMAX/50.) • FLOAT(I - 1) + VO 101 30 CONTINUE C c VMAX = VMAX • UTH C C AXES LENGTHS FOR THE PLOTS C AXL1 =5.0 AXL2 = 7.0 AXL3 = 5.0 C C INIT I A L I S E VELOCITY RANGES IN PHASE SPACE C PLOT. AS WELL THE PARTICLE SAMPLING INTERVAL C I SAMP IS DETERMINED FROM THE TOTAL PARTICLE C NUMBER AND THE SAMPLE SIZE NSAMP. C IF (VMAX2.EQ.0.) VMAX2=VMAX/UTH VMINPH=-VMAX2 VMAXPH=VMAX2 ISAMP=NP/NSAMP IF (NP.LT.NSAMP) ISAMP=1 C C C C C C C C MSHSEL RECORDS WHICH MESH POINTS YOU WANT TO USE C IN THE TIME HISTORY OF THE FIELDS. NMSEL IS THE C NUMBER OF MESH POINTS YOU'RE INTERESTED IN. C IF NMSEL=0 NO MESH POINTS ARE WANTED C IF NMSEL = -1 EVERY MESH POINT IS WANTED C SOMETHING SIMILAR IS DONE FOR NFTSEL AND MFTSEL C C IF (INPSW(11) .EQ. 1) GO TO 50 IF (NMSEL .GE. 0) GO TO 50 C DO 40 I = 1, MESH MSHSEL(I) = I 40 CONTINUE C NMSEL = MESH 50 CONTINUE C IF (NFTSEL .GE. 0) GO TO 70 C DO 60 I = 1, M22 MFTSEL(I) = I 60 CONTINUE C NFTSEL = M22 102 C C C-C C C C C c c c c c c-c c c c 70 80 90 130 CONTINUE IF (NUMTR .EQ. 0) GO TO 80 WRITE (3,170) NUMTR IF (NFTSEL .EQ. 0) GO TO 90 WRITE (17,130) INVLFT WRITE (17,140) WRITE (17,150) (MFTSEL(I),1=1,NFTSEL) RETURN FORMAT ('DATA CONSISTS OF THE FOURIER TRANSFORM VALUES' 'OF THE EX FIELD FOR EVERY ', 13, ' TIME STEPS') ('A'f 'THE K NUMBERS RECORDED ARE ') (2013) ('TRACER POSITIONS IN TERMS OF PLASMA LENGTH'/, 140 FORMAT 150 FORMAT 170 FORMAT 1 2 END 'AND TRACER ' TRACERS') VELOCITIES, IN VTH UNITS. THERE '/, SUBROUTINE PREDGN(NTS, NTSC) IMPLICIT REAL•4(A - H,0 - Z,X) COMMON /PARAR/ XP(100000), UX(100000) COMMON /FLDAR/ EX( 1000), RHOdOOO), RHODOOOO) COMMON /EMCON/ EOVERM, QET, QDX, C, PI COMMON /DENCON/ OMEGAP, UTH, CID, DENSUR COMMON /LASCON/ OMEGA, ALAMBA, TIME COMMON /DGNS/ VMAX, INPSW(30), MSHSEL(200), NMSEL, 1 MFTSEL(129), NFTSEL, INVLPL, INVLFT, 2 INVLFL, INVLEN, INVLXV COMMON /PARC/ A1, ANE1, PHI 1 , A2, ANE2, PHI 2, NTPULS, 1 NUMTR, TRP(30), TRV(30), INVLTR, UTHIN, OMEGIN, 2 ANL, VO COMMON /MSHCON/ MESH, MESHM1, NP, DX, XL, DT TS = FLOAT(NTS + NTSC) TIMP IS THE TIME IN PLASMA PERIODS. TIMP = TS • DT • OMEGAP / (2.«PI) 11 CONTINUE IF (MOD(NTS,INVLPL).EQ.0) CALL MYPLOT(UTH, VO, TIMP) INVLFL OUTPUT=16 INVLFT OUTPUT=17 INVLEN 0UTPUT=15 PHASPC PLOTS A PHASE SPACE SCATTER PLOT OF THE ELECTRONS. 103 IF (MOD(NTS,INVLXV).EQ.O) .CALL PHASPC(UTH,XL,NP,TIMP,INPSW(10)) C C C c C TRACER HISTORY SECTION C c c c-IF (NUMTR .EQ. 0) GO TO 20 IF (MOD(NTS,INVLTR) .EQ. 0) .CALL TRACER(NUMTR, UTH) 20 CONTINUE IF (NMSEL .EQ. 0) GO TO 30 C IF (MOD(NTS,INVLFL) .EQ. 0) CALL RMSF(MESH) 30 CONTINUE C c c c C FOURIER TRANSFORM HISTORY SECTION C IF (NFTSEL .EQ. 0) GO TO 40 IF (MOD(NTS,INVLFT) .EQ. 0) CALL FTRANS(MESH) 40 CONTINUE C c c c C ENERGY CALCULATES TOTAL ENERGY C TOTAL ENERGY=KINETIC ENERGY+ELECTROSTATIC C ENERGY C PLACES K.E., EL.E. A TOTAL E. IN DEVICE 15 IF (MOD(NTS,INVLEN) .EQ. 0) .CALL ENERGY(DENSUR,UTH) C C c-RETURN END C c 0========================================================= c C ENERGY CALCULATES THE ENERGY IN THE SYSTEM C SUBROUTINE ENERGY(DENSUR,UTH) COMMON /STRGE/ VOLD(100000),XOLD(100000),RHOLD(1000) COMMON /FLDAR/ EX(lOOO), RHO(1000), RHOD(lOOO) 104 COMMON /MSHCON/ MESH, MESHM1, NP, DX, XL, DT COMMON /DGNS/ VMAX, INPSW(30), MSHSEL(200), NMSEL, 1 MFTSEL(129), NFTSEL, INVLPL, INVLFT, 2 INVLFL, INVLEN, INVLXV C C THE ELECTROSTATIC ENERGY IS FOUND BY INTEGRATING C OVER X. THE METHOD USED IS SIMPSONS C AS SIMPSONS METHOD REQUIRES AN EVEN NUMBER OF PANELS C AND WE ALWAYS HAVE AN ODD NUMBER, THE LAST PANEL IS C ADDED ON SEPARATELY (ENDVAL) DX3 = DX / 3. EN = 0. C DO 10 I = 3, MESHM1, 2 SQ = EX(I - 2).«2 + 4. -(EX(I - 1 ) " 2 ) + EX(l)-«2 EN = EN + SQ • DX3 10 CONTINUE C ENDVAL = (EX(MESHM1)••2 + EX(MESH).-2) • .5 • DX ENEL = (EN + ENDVAL) / 25.1327 EN = 0. C C 25.1327=8-PI C DO 20 I = 1, NP EN = EN + VOLD(I) •• 2 20 CONTINUE C C C 4.5535E-28 = 0.5 • MASS OF ELECTRON(CGS UNITS) C VELRMS = SQRT(EN/FLOAT(NP))/UTH ENKIN = EN • 4.5535E-28 • DENSUR ENTOT = ENEL + ENKIN WRITE (15) ENKIN, ENEL, ENTOT, VELRMS RETURN END C C . c c C XAVER INTEGRATES THE ARRAY DIN OVER X VIA C SIMPSONS RULE.V C SUBROUTINE XAVER(DIN, X) DIMENSION DIN(1) COMMON /MSHCON/ MESH, MESHM1, NP, DX, XL, DT XTO = 0 . C C C SIMPSON'S INTEGRATION METHOD IS USED HERE. BUT C SIMPON'S METHOD REQUIRES AN EVEN NUMBER OF PANELS 105 C OR CELLS. IN THIS PROGRAM MESH IS ALWAYS EVEN AND C HENCE THERE ARE 'AN ODD NUMBER OF CELLS. TO GET C AROUND THIS, THE LAST CELL IS NOT INCLUDED IN C THE INTEGRATION AND ITS AREA IS ADDED ON AFTER. C DX3 = DX / 3. DO 10 I = 3, MESHM1, 2 XTO = XTO+(DIN(l - 2) + 4.«DIN(I - 1)+DIN(I))-DX3 10 CONTINUE ENDVAL = (DIN(MESHM1) + DIN(MESH)) • .5 • DX X = (XTO + ENDVAL) / XL RETURN END C C TRACER WRITES INTO DEVICE 3 THE PARTICLE C POSITION AND VELOCITY IN TERMS OF XL AND UTH. SUBROUTINE TRACER(NUMTR, UTH) COMMON /STRGE/ VOLD(100000),XOLD(100000),RHOLD(1000) COMMON /MSHCON/ MESH, MESHM1, NP, DX, XL, DT C C NPP = NP - NUMTR + 1 C DO 10 I = NPP, NP XOLD(I) = XOLD(I) / X L VOLD(I) = VOLD(I) / UTH 10 CONTINUE WRITE (3) (XOLD(I),I=NPP,NP), (VOLD(I),I=NPP,NP) DO 20 I = NPP, NP XOLD(I) = XOLD(I) • XL VOLD(I) - VOLD(I) • UTH 20 CONTINUE C RETURN END C C c C RMSF CALCULATES THE RMS E FIELD AND THE RMS CHARGE C DENSITY C IT CAN ALSO STORE THE VALUES OF THESE TWO ARRAYS C FOR DIFFERENT MESH POINTS C C SUBROUTINE RMSF(MESH) REAL* 4 STOR1O 000), STOR2(1000) COMMON /STRGE/ VOLD(100000),XOLD(100000),RHOLD(1000) C C 106 C c c c c COMMON /FLDAR/ EX(lOOO), RHOO 000), RHODOOOO) COMMON /DGNS/ VMAX, INPSW(30), MSHSEL(200), NMSEL, 1 MFTSEL(129), NFTSEL, INVLPL, INVLFT, 2 INVLFL, INVLEN, INVLXV NSTOR = 0 10 IF (INPSW(11) .EQ. 0) GO TO 30 DO 20 I = 1, MESH STOR1(I) = E X ( I ) •• 2 STOR2(l) = RHOLD(I) •• 2 20 CONTINUE CALL XAVER(STOR1, AV) AV = SQRT(AV) STOR1(1) = AV CALL XAVER(STOR2,.AV) STOR1(2) = SQRT(AV) NSTOR = 2 IF (INPSW(11) .EQ. 1) GO TO 50 30 CONTINUE DO 40 I. = 1 , NMSEL I I = I + NSTOR S T G R K I I ) = EX(MSHSEL(I ) ) NST = I I + NMSEL STORI(NST) = RHOLD(MSHSEL(I)) 40 CONTINUE 50 CONTINUE WRITE (16) (STOR1(I),1=1,NST) RETURN END C C c C FTRANS FOURIER TRANSFORMS EX (AND RHO) C AND WRITES INTO 17 THE VALUES OBTAINED OF C THE MODES SELECTED BY MFTSEL. SUBROUTINE FTRANS(MESH) REAL•4 FT(512) COMMON /STRGE/ VOLD(100000),XOLD(100000),RHOLD(1000) COMMON /FLDAR/ EX(lOOO), RHO(lOOO), RHOD(1000) COMMON /DGNS/ VMAX, INPSW(30), MSHSEL(200), NMSEL, 1 MFTSEL(129), NFTSEL, INVLPL, INVLFT, 2 INVLFL, INVLEN, INVLXV C C WARNING*** —>WHEN CHOOSING A PARTICULAR 107 C K NUMBER=I, TO TAKE ITS HISTORY YOU PICK C FT(I+1) AS THERE IS NO FT(0) IN FORTRAN. C K=0 CORRESPONDS TO FT(1) C THAT'S WHY THERE IS MFTSEL(l)+1 C M2 = MESH / 2 CALL FTRXK(EX, FT, M2) WRITE (17) (FT(MFTSEL(I) + 1),I=1,NFTSEL) C CALL FTRXK(RHOLD,FT,M2) C WRITE (15) (FT(MFTSEL(I)+1),1=1,NFTSEL) RETURN END C c c C MYPLOT PLOTS EX IT FOURIER TRANSFORM, RHO, ITS C FOURIER TRANSFORM AND THE VELOCITY DISTRIBUTION. C C SUBROUTINE MYPLOT(UTH,VO, TIMP) COMMON /PLCOM/ X I ( 1 0 0 0 ) r X2(1000), X 3 O 0 0 0 ) , AXL1, AXL2, AXL3 COMMON /MSHCON/ MESH, MESHM1, NP, DX, XL, DT COMMON /FLDAR/ EX(1000), RHO( 1000), RHODOOOO) COMMON /STRGE/ VOLD(100000),XOLD(100000),RHOLD(1000) COMMON /DGNS/ VMAX, INPSW(30), MSHSEL(200), NMSEL, . 1 MFTSEL(129), NFTSEL, INVLPL, INVLFT, 2 INVLFL, INVLEN, INVLXV INTEGER.4 TLE(3) /'VEL.', ' EX *, 'RHO '/ INTEGER.4 XAXT(3) /'VEL.', ' X ', ' K '/ INTEGER* 4 YAXT(3) / ' F ( V ) ' , 'EX ', 'RHO '/ REAL.4 F T ( 6 0 0 ) , V(101) M2 = MESH / 2 IPLL = 0 C 10 IPLL = IPLL + 1 CALL MYZERO(FT, M2) IF ( I P L L .EQ. 4) GO TO 90 IF (INPSW(IPLL) .EQ. 0) GO TO 10 IF (INPSW(IPLL) .EQ. 2) GO TO 20 C C HERE THE USER IS PROMPTED IF HE WANTS A PLOT. C WRITE (6,100) T L E ( I P L L ) READ (5,110) ISWPL IF (ISWPL .EQ. 0) GO TO 10 20 GO TO (30, 70, 8 0 ) , IPLL C C CALCULATE THE VELOCITY DISTRIBUTION C DONE IN THE USUAL MANNER, WITH BINS ( V ( J ) ) . 30 MESHV = 51 MESHV1 = MESHV - 1 C 108 C c c CALL MYZERO(V, MESHV) DV = VMAX / FLOAT(MESHV1) DO 50 I = 1, NP DO 40 J = 1, MESHV XV = DV • FLOAT(J) UT = ABS(VO-UTH-ABS(VOLD(I))) IF (UT .GE. XV) GO TO 40 V ( J ) = V ( J ) + 1. C J = MESHV GO TO 45 40 CONTINUE C 45 CONTINUE 50 CONTINUE C c DVH = VMAX / (UTH*FLOAT(MESHV1)) C C TO NORMALIZE DISTRIBUTION ( F ( V ) ) , V ( j ) IS C DIVIDED BY THE NUMBER OF OCCURENCES (NP), THE C BIN WIDTH (DVH), AND 2 SINCE INSTEAD OF PLACING C VEL. IN A BIN FROM -VMAX TO +VMAX, ABS(VEL.) WAS C PLACED IN A BIN FROM 0 TO +VMAX C DO 60 J = 1 , MESHV V ( J ) = V ( J ) / (FLOAT(NP)•DVH•2.) 60 CONTINUE C C c c c c c c c c c c c c c c c-c c CALL GPLOT(V, X3, MESHV, XAXT(1), YAXT(1), AXL3, TIMP) GO TO 10 THIS BIT CHECKED OUT IF F(V) WAS NORMALIZED (AREA=1) WRITE(4,340) (V(J),J=1,MESHV) XTI = 0. DO 111 I = 2, MESHV1,2 XTI = XTI + (V(I - 1) + 4 . - V ( l ) + V(I+1))-DVH/3. 111 CONTINUE XTI=2.-XTI WRITE (6,121) XTI 121 FORMAT(1X,'AREA OF DISTIBUTION=',F9 .4) c 109 C c c c 70 CALL GPLOT(EX, X1 , MESH, XAXT(2), YAXT(2), AXL1, TIMP) C CALL FTRXK(EX, FT, M2) C CALL GPLOT(FT, X2, M2, XAXT(3), YAXT(2), AXL2, TIMP) GO TO 10 C C 80 CALL GPLOT(RHOD, X1, MESH, XAXT(2), YAXT(3), AXL1, TIMP CALL FTRXK(RHOD, FT, M2) CALL GPLOT(FT, X2, M2, XAXT(3), YAXT(3), AXL2, TIMP) GO TO 10 C c 90 CONTINUE RETURN 100 FORMAT (/'DO YOU WANT ', A4, ' PLOT7 (0=NO 1=YES)') 110 FORMAT (11) END C C C c C SIMPLE PLOTTING ROUTINE C SUBROUTINE GPLOT(RAV, RAH, NDIM, XAXT, YAXT, AHL, TIMP) INTEGER-4 YAXT, XAXT DIMENSION RAV(NDIM), RAH(NDIM), DHIN(550), DVIN(550) C C DO 10 I I = 1, NDIM DHIN(II) = RAH(II) D V I N ( I I ) = R A V ( I I ) 10 CONTINUE C VL = 5. CALL SCALE(DVIN, NDIM, VL, YMIN, DY, 1) CALL SCALE(DHIN, NDIM, AHL, XMIN, DX, 1) CALL AXIS(0., 1., XAXT, -4, AHL, 0., XMIN, DX) CALL AXIS(0., 1., YAXT, 4, VL, 90., YMIN, DY) C CALL LINE(DHIN,DVIN,NDIM,1) CALL PLOT( (DHIN( 1 ) - 0 . ) / L , (DVIN(l) + 1.)/1., 3) C DO 20 J J = 2, NDIM CALL PLOT( (DHIN(JJ) - 0 . ) / L , (DVIN(JJ) + 1.)/1., 2) 20 CONTINUE C CALL NUMBER(2., 7., .14, TIMP, 0., 3) CALL PLOT(0.00, 6.00, 3) CALL PLOT(0.00, 10.0, 2) CALL PLOT(7.0, 10.0, 2) 110 CALL P L O T ( 7 . 0 , 1 . 0 , 2 ) CALL PLOT(AHL, 1.0, 2) CALL PLOTO0.5, 10.5, 3) CALL PLOTO0.5, 10.0, 2) CALL PLOT(15., 0., -3) RETURN END C C C c c SUBROUTINE FTRXK(RIN, FT, MOH) DIMENSION YA(600), Y S ( 6 0 0 ) , R I N ( 1 ) , FT(1) C C THIS SUBROUTINE CALCULATES THE FOURIER C TRANSFORM OF RIN (WITHOUT CHANGING RIN) C AND PLACES IT IN THE ARRAY FT C THE DIMENSION OF RIN IS MESH C THE DIMENSION OF FT IS MOH C I F MESH=2•-N (N BEING AN INTEGER FROM 5 TO C 9 MESH=32,64,128,256,512) THE FFT C ROUTINES ( I . E . C64,S64) ARE USED. IF THIS C IS NOT THE CASE FOUR2 IS USED C IN BOTH CASES THE OUTPUT FT IS THE MODULUS C SQRT(FT X FT-) MOHM1 = MOH - 1 MOHP1 = MOH + 1 MESH = MOH • 2 IF (MOH .LE. 256) GO TO 40 C 10 CONTINUE C DO 20 I = 1, MESH F T ( I ) = R I N ( I ) 20 CONTINUE C CALL FOUR2(FT, MESH, 1, 1, 0) C DO 30 I = 1, MOHM1 J = 2 • I + 1 TS = F T ( J ) •• 2 TA = F T ( J + 1) •• 2 FT(I+1) = SQRT(TS + TA) / FLOAT(MOH) 30 CONTINUE FT(1)=0. C GO TO 110 C C 40 CONTINUE C C FOR THE FFT ROUTINES THE DATA IS PREPARED INTO 111 C A SYMMETRIC ARRAY (YS) AND A ANTISYMMETRIC ARRAY C (YA) SEE UBC WRITEUP FFT YA(1) = 0 . YS(1) = RIN(1) YS(M0HP1) = RIN(MOHPI) C DO 50 J = 1, MOHM1 NP = MOHP1 + J NM = MOHP1 - J YA(NM) = -.5 • (RIN(NP) - RIN(NM)) YS(NM) = .5 • (RIN(NP) + RIN(NM)) 50 CONTINUE C C IF (MOH .NE. 32) GO TO 60 CALL C64(YS, YS) CALL S64(YA, YA) GO TO 90 C C 60 IF (MOH .NE. 64) GO TO 70 CALL C128(YS, YS) CALL S128(YA, YA) GO TO 90 C C 70 IF (MOH .NE. 128) GO TO 80 CALL C256(YS, YS) CALL S256(YA, YA) GO TO 90 C C 80 IF (MOH .NE. 256) GO TO 10 CALL C512(YS, YS) CALL S512(YA, YA) C C 90 CONTINUE C C DO 100 J = 2, MOH TA = YA(J) • • 2 TS = YS(J) •• 2 F T ( J ) = SQRT(TA + TS) / FLOAT(MOH) 100 CONTINUE C FT(1) = 0. C C 110 CONTINUE RETURN END 112 C C c c SUBROUTINE PHASPC(VTH,XL,NP,TIMP,ISTRM) C COMMON /STRGE/ VOLD(100000),XOLD(100000),RHOLD(1000) COMMON /PHASPL/ VMAX2,VMIN,VMAX,I SAMP,NSAMP C C REAL•4 OFFH,OFFL,VS(100000),XS(100000) C AXH=10. AXV=6. XO=0.0 YO=0.0 TOL=.1 ITER=0 C C CHECK I F ITS THE STREAMING CASE ISTRM=3 IF (ISTRM.NE.3) GO TO 11 C C TO SAVE PLOTTING TIME, THE SAMPLE POINTS ARE C SORTED SO THAT ONE BEAM IS PLOTTED FIRST C AND THEN THE OTHER. J=0 DO 12 I = 1,NP,I SAMP IF (VOLD(l).LT.0.) GO TO 12 J=J+1 VS(J)=V0LD(I)/VTH XS(J)=XOLD(l)/XL 12 CONTINUE DO 13 1=1,NP,ISAMP IF (VOLD(I).GE.0.) GO TO 13 J=J+1 VS(J)=VOLD(l)/VTH XS(J)=X0LD(I)/XL 13 CONTINUE NPHAS=J GO TO 14 11 CONTINUE J=0 DO 15 1=1,NP,ISAMP J=J+1 VS(J)=VOLD(l)/VTH XS(J)=XOLD(l)/XL 15 CONTINUE NPHAS=J C 14 CONTINUE 113 5 OFFL=0. OFFH=0. AVH=0. AVL=0. C DO 10 1=1,NPHAS C IF (VS(I).LE.VMAX) GO TO 19 OFFH=OFFH+1. AVH=AVH+VS(I) GO TO 10 19 LF (VS(I).GE.VMIN) GO TO 10 OFFL=OFFL+1. AVL=AVL+VS(I) 10 CONTINUE C C IF MORE THAN TOL OF THE POINTS ARE GOING OUTSIDE C THE VELOCITY RANGES, VMAX ON THE PLOT IS INCREASED. C TO A VALUE SLIGHTLY LARGER THAN THE VELOCITY OF THE C FASTEST PARTICLE. C PERH=OFFH/FLOAT(NP) PERL=OFFL/FLOAT(NP) IF ((PERH.LT.TOL).AND.(PERL.LT.TOL)) GO TO 20 IF (PERH.LT.TOL) GO TO 30 AVH=AVH/OFFH VMAX=FLOAT(INT(AVH•1.414)) +1. 30 AVL=AVL/OFF L VMIN=FLOAT(lNT(AVL.1.414)) - 1 . ITER=ITER+1 IF (ITER.LT.3) GO TO 5 C C 20 DX=1./AXH VWID=ABS(VMAX-VMIN) DV=VWID/AXV 'CALL AXIS(XO,YO,*L' ,-1 ,AXH,0.0,0.0,DX) CALL AXIS(XO,YO,'VTH',3,AXV,90.,VMIN,DV) C DO 40 1=1,NPHAS IF ((VS(I).GT.VMAX).OR.(VS(I).LT.VMIN)) GO TO 40 XPL=XS(I)•AXH + XO VPL= (ABS(VS(I)-VMIN)/VWID)•AXV + YO CALL SYMBOL(XPL,VPL,.07,14,0.0,-1) 40 CONTINUE C PERL=PERL•100. PERH=PERH-100. XN=AXH-1.+XO YNL= YO+.1 YNH=YO+AXV-.3 CALL NUMBER(XN,YNL,.14,PERL,0.0,1) CALL NUMBER(XN,YNH,.14,PERH,0.0,1) 1 14 CALL NUMBER(XN,YNH+.4, .14, TIMP, 0., 3) CALL PLOT(XO,AXV+YO,3) CALL PLOT(XO+AXH,AXV+YO,2) CALL.PLOT(XO+AXH,YO,2) CALL SYMBOL(0.,0.,.02,3,0.,~1) CALL SYMBOL(12.7,9.7,.02,3,0.,-1) CALL PLOT(15.,0.,-3) C c RETURN END 115 PARTICLE MOVING ROUTINE ASSEMBLER VERSION OF MOVEP (STRICTLY ELECTROSTATIC) FRO EQU FR2 EQU FR4 EQU FR6 EQU EXTRN EDGE ENTRY MOVEPA MOVEPA ENTER 12,SA=SAVE1 INITIALIZATION • L 2,=A(XP) L 3,=A(UX) L 4,=A(EX) L 5,=A(RHO) L 11,=A(MESHM1) L 7,0(11) M 6,=F'4' AR 7,5 LR 6,7 S 6,=F'4* L 10,=A(DT) LE FRO,0(10) STE FRO,DTT L 10,=A(DX) LE FRO,0(10) STE FR0,DXX L 10,=A(XL) LE FR0,0(10) STE FRO,XLL L 10,=A(QET) LE FRO,0(10) STE FR0,QETT L 10,=A(UTH) LE FRO,0(10) STE FRO,UTHH LE FR0,=E'1.0' STE FRO,PSIGN LE FRO,=E'-1.0' STE FRO,MSIGN L 10,=A(CID) L 9,=A(QDX) LE FRO,0(9) STE FR0,QDXX ME FR0,0(10) STE FRO,CION L 10,=A(NP) L 9,0(10) M 8,=F'4' S 9,=F'4' 7=ADDRESS OF 6=ADDRESS OF GET SOME FROM THE RHO(MESH) RHO(MESHM1) OF THE CONSTANTS COMMON BLOCKS AND PLACE INTO LOCAL STORAGE. 9=NP 1 16 AR 9,2 LE FRO,=E'0.0' STE FRO,RHE PARTICLE LOOP LOOP1 LE FR4,0(2) DE FR4,DXX EFIX FR4,11 « STD FR4,INT0 THIS IS THE RETURN POINT FOR THE MAIN LOOP. THE NEXT FEW LINES ARE THE EQUIVALENT OF J=INT(X/DX) NOTE THAT 1 IS NOT ADDED TO J AS IS DONE IN THE FORTRAN VERSION SINCE THE FIRST ARRAY ELEMENT IN ASSEMBLER IS REFERRED TO BY THE ADDRESS OF THE ARRAY. LD FR4,=D'0' LE FR4,INT0 AW FR4,=X'4E00000080000000' STD FR4,INT0+8 XI INT0+12,X'80' L 11,INT0+12 LD FR4,INT0 FLOAT (11),FR6 THE NEXT GROUP OF LINES R=FLOAT(J) ST 11,FLC+4 XI FLC+4,X'80' LD FR6,FLC SD FR6,=X'4E00000080000000' FR4=PX=XTEST-R FR6=PX1 EX(J)-PX1 EX(J1).PX FR4=EXI SER FR4,FR6 LCER FR6,FR4 AE FR6,=E'1.0' M 10,=F'4' AR 11,4 ME FR6,0(11) ME FR4,4(11) AER FR4,FR6 ME FR4,QETT AE FR4,0(3) STE FR4,0(3) ME FR4,DTT AE FR4,0(2) CE FR4,XLL BNH JUMP1 STE FR4,XPT LE FR4,0(3) STE FR4,UXT CALL EDGE,(XPT,UXT,XLL,UTHH,MSIGN), VL LE FR4,UXT STE FR4,0(3) LE FR4,XPT FIND THE ELECTRIC FIELD AT THE PARTICLE POSITION. DETERMINE THE NEW VELOCITY DETERMINE THE NEW POSITON. SEE IF PARTICLE HAS GONE PAST RIGHT EDGE. IF SO CALL FORTRAN ROUTINE EDGE 1 17 JUMP1 .0' B JUMP2 CE FR4,=E'0 BNL JUMP2 STE FR4,XPT LE FR4,0(3) STE FR4,UXT CALL EDGE,(XPT,UXT,XLL,UTHH,PSIGN),VL LE FR4,UXT STE FR4,0(3) LE FR4,XPT JUMP2 STE FR4,0(2) DE FR4,DXX CHECK TO SEE IF PARTICLE HAS PAST LEFT EDGE. IF SO CALL EDGE. EFIX FR4,11 STD FR4,INT0 LD FR4,=D'0' LE FR4,INT0 AW FR4,=X'4E00000080000000' STD FR4,INT0+8 XI INT0+12,X'80' L 11,INT0+12 LD FR4,INT0 FLOAT (11),FR6 ST 11,FLC+4 • XI FLC+4,X'80' LD FR6,FLC SD FR6, =X'4E00000080000000' THIS SECTION CALCULATES THE CHARGE DENSITIES (RHO) AT THE MESH POINTS. FIRST DETERMINE WHERE PARTICLE I S . PX AND PX1 DETERMINED JUST AS THEY WERE WHEN FINDING EXI SER FR4,FR6 LCER FR6,FR4 AE FR6,=E'1.0' M 10,=F'4' AR 11,5 ME FR4,QDXX ME FR6,QDXX CR 11,7 BE JUMP3 CR 11,5 BNE JUMP4 AE FR4,4(11) STE FR4,4( 1 1 ) JUMP3 AE FR6,RHE STE FR6,RHE B JUMP6 JUMP4 CR 11,6 BNE JUMP5 AE FR4,RHE STE FR4,RHE AE FR6,0(11) STE FR6,0(11) FR4=PX FR6=PX1 FR4=PX-QDX FR6=PX1.QDX J=MESH7 IF J=1,MESHM1 OR MESH CHARGE CONTRIBUTION GOES TO RHE. J=1 7 J=MESHM17 B JUMP6 JUMP5 AE FR4,4(11) STE FR4,4(11) AE FR6,0(1 1 ) STE FR6,0(11) JUMP6 LA 2,4(2) LA 3,4(3) CR 2,9 BNH L00P1 PATH COMES 1<J<MESHM1 HERE IF INCREMENT PARTICLE ADDRESS. I=NP7 OR K N P GO BACK TO LOOP1 END OF MAIN PARTICLE LOOP CALCULATE RHO LE FR4,RHE STE FR4,0(5) STE FR4,0(7) LOOP2 LE FR4,0(5) SE FR4,CION STE FR4,0(5) LA 5,4(5) CR 5,7 BNH LOOP2 SUBTRACT CONSTANT ION CHARGE (CION) FROM RHO AT EACH MESH POINT EXIT LTORG SAVE 1 DS 18F PSIGN DS E MSIGN DS E UXT DS E XPT DS E UTHH DS E CION DS E RHE DS E QDXX DS E QETT DS E DTT DS E XLL DS E DXX DS E INTO DC 2D'0' FLC • DC X ' 4 E 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ' • FLDAR COM EX DS 1000E RHO DS 1 000E RHOD DS 1000E PARAR COM XP DS 100000E UX DS 100000E • MSHCON COM MESH DS F MESHM1 DS F NP DS F DX DS E XL DS E DT DS E • EMCON COM EOVERM DS E QET DS E QDX DS E C DS E PI DS E • DENCON COM OMEGAP DS E UTH DS E CID DS E DENSUR DS E END 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            data-media="{[{embed.selectedMedia}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0085213/manifest

Comment

Related Items