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