1 -- 2 -- Patricia Fasel 3 -- Los Alamos National Laboratory 4 -- 1990 August 5 -- 6 -- Copyright, 1990, The Regents of the University of California. 7 -- This software was produced under a U.S. Government contract (W-7405-ENG-36) 8 -- by the Los Alamos National Laboratory, which is operated by the University 9 -- of California for the U.S. Department of Energy. The U.S. Government is 10 -- licensed to use, reproduce, and distribute this software. Permission is 11 -- granted to the public to copy and use this software without charge, provided 12 -- that this notice and any statement of authorship are reproduced on all 13 -- copies. Neither the Government nor the University makes any warranty, 14 -- express or implied, or assumes any liability for the use of this software. 15 16 module Pic (pic) where 17 18 import PicType 19 import Consts 20 import Utils 21 import ChargeDensity 22 import Potential 23 import ElecField 24 import PushParticle 25 import Array--1.3 26 27 -- PIC, particle in cell, a basic electrodynamics application 28 -- Given an initial configuration of particles, follow how they move under the 29 -- electric field they induce 30 -- Torroidal boundary conditions are assumed, so wrap in both dimensions 31 -- given nPart the number of particles considered 32 -- given nStep the number of time steps to put the particles through 33 -- given nCell the dimension of the matrix of cells 34 35 pic :: Indx -> [Char] 36 pic nPart = 37 show dt' 38 where 39 partHeap = initParticles nPart 40 dt = 0.001 41 phi = initPhi partHeap 42 (dt', phi', partHeap') = timeStep partHeap phi dt 0 nStep 43 44 45 -- during each time step perform the following calculations 46 -- calculate the charge density (rho), using position of particles 47 -- calculate the new potential (phi), by solving Laplace's equation 48 -- del2(phi) = rho , using rho and phi of last timestep 49 -- calculate the electric field, E = del(phi), using phi of this time step 50 -- push each particle some distance and velocity using electric field, for a 51 -- timestep deltaTime, small enough that no particle moves more than the 52 -- width of a cell 53 -- an NxN mesh is used to represent value of x and y in the interval [0,1] 54 -- so delta_x = delta_y = 1/n 55 -- 56 -- phi ((0,0), (n,n)) = electrostatic potential at grid point (i,j) 57 -- rho ((0,0), (n,n)) = charge density at grid point (i,j) 58 -- xElec ((0,0), (n,n)) = x component of electric field between (i,j) (i,j+1) 59 -- yElec ((0,0), (n,n)) = y component of electric field between (i,j) (i+1,j) 60 -- [xyPos] = (x,y) coordinate of particle displacement in units of delta_x 61 -- [xyVel] = (x,y) coordinate of particle velocity in units of delta_x/sec 62 63 timeStep :: ParticleHeap -> Phi -> Value -> Indx -> Indx -> (Value, Phi, ParticleHeap) 64 timeStep partHeap phi dt depth step 65 | step == 0 = (dt, phi, partHeap) 66 | otherwise = 67 timeStep partHeap' phi' dt' depth' (step-1) 68 where 69 rho = chargeDensity partHeap 70 phi' = potential phi rho depth 1 71 xyElec = elecField phi' 72 (maxVel, maxAcc, partHeap') =pushParticle partHeap xyElec dt 0 0 73 dt' = (sqrt (maxVel*maxVel + 2*maxAcc) - maxVel) / maxAcc 74 depth' = (depth+1) `rem` maxDepth 75 76 77 initParticles :: Indx -> ParticleHeap 78 initParticles nPart = 79 (xyPos, xyVel) 80 where 81 nCellD = fromIntegral nCell 82 nPartD = fromIntegral (nPart+1) 83 xyPos = [(xPos,yPos) | 84 i <- [1..nPart], 85 xPos <- [nCellD * genRand (fromIntegral i/ nPartD)], 86 yPos <- [nCellD * genRand xPos]] 87 xyVel = [(0.0,0.0) | i <- [1..nPart]] 88 89 90 initPhi :: ParticleHeap -> Phi 91 initPhi partHeap = 92 potential phi0 rho maxDepth 1 93 where 94 rho = chargeDensity partHeap 95 phi0 = array ((0,0), (n,n)) 96 [((i,j), 0.0) | i <- [0..n], j <- [0..n]] 97 n = nCell-1