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