1 -- 
    2 --      Patricia Fasel
    3 --      Los Alamos National Laboratory
    4 --      1990 August
    5 --
    6 module ChargeDensity (chargeDensity) where
    7 
    8 import  PicType
    9 import  Consts
   10 import Array--1.3
   11 
   12 -- Phase I: calculate the charge density rho
   13 -- Each particle represents an amount of charge distributed over an entire cell.
   14 -- In discretizing the charge density at the grid points we split the cell into
   15 --    four rectangles and assign to each corner an amount of charge proportional
   16 --    to the area of the opposite diagonal sub rectangle
   17 -- So each particle contributes a portion of charge to four quadrants
   18 -- particle position is (i+dx, j+dy)
   19 -- nw quadrant (1-dx)(1-dy)(charge) is added to rho(i,j)
   20 -- ne quadrant (dx)(1-dy)(charge) is added to rho(i,j+1)
   21 -- sw quadrant (1-dx)(dy)(charge) is added to rho(i+1,j)
   22 -- se quadrant (dx)(dy)(charge) is added to rho(i+1,j+1)
   23 -- wrap around used for edges and corners
   24 
   25 chargeDensity :: ParticleHeap -> Rho
   26 chargeDensity (xyPos, xyVel) =
   27         accumArray (+) 0 ((0,0), (n,n)) (accumCharge xyPos)
   28         where
   29             n = nCell-1
   30 
   31 
   32 -- for each particle, calculate the distribution of density
   33 -- based on (x,y), a proportion of the charge goes to each rho
   34 
   35 accumCharge :: [Position] -> [MeshAssoc]
   36 accumCharge [] = []
   37 accumCharge ((x,y):xys) =
   38         [((i ,j ) , charge * (1-dx) * (1-dy))] ++
   39         [((i',j ) , charge * dx * (1-dy))] ++
   40         [((i ,j') , charge * (1-dx) * dy)] ++
   41         [((i',j') , charge * dx * dy)] ++
   42         accumCharge xys
   43         where
   44             i = truncate x
   45             i' = (i+1) `rem` nCell
   46             j = truncate y
   47             j' = (j+1) `rem` nCell
   48             dx = x - fromIntegral i
   49             dy = y - fromIntegral j