{-# OPTIONS -Wall #-}

{- | 
Module      :  LPFPCore.Maxwell
Copyright   :  (c) Scott N. Walck 2023
License     :  BSD3 (see LICENSE)
Maintainer  :  Scott N. Walck <walck@lvc.edu>
Stability   :  stable

Code from chapter 29 of the book Learn Physics with Functional Programming
-}

module LPFPCore.Maxwell where

import LPFPCore.SimpleVec
    ( R, Vec(..), (^/), (^+^), (^-^), (*^)
    , vec, negateV, magnitude, xComp, yComp, zComp, iHat, jHat, kHat )
import LPFPCore.CoordinateSystems
    ( ScalarField, VectorField
    , cart, shiftPosition, rVF )
import LPFPCore.ElectricField ( cSI, mu0 )
import qualified Data.Map.Strict as M

directionalDerivative :: Vec -> ScalarField -> ScalarField
directionalDerivative :: Vec -> ScalarField -> ScalarField
directionalDerivative Vec
d ScalarField
f Position
r
    = (ScalarField
f (Vec -> Position -> Position
shiftPosition (Vec
d Vec -> R -> Vec
^/ R
2) Position
r) R -> R -> R
forall a. Num a => a -> a -> a
- ScalarField
f (Vec -> Position -> Position
shiftPosition (Vec -> Vec
negateV Vec
d Vec -> R -> Vec
^/ R
2) Position
r))
      R -> R -> R
forall a. Fractional a => a -> a -> a
/ Vec -> R
magnitude Vec
d

curl :: R -> VectorField -> VectorField
curl :: R -> VectorField -> VectorField
curl R
a VectorField
vf Position
r
    = let vx :: ScalarField
vx = Vec -> R
xComp (Vec -> R) -> VectorField -> ScalarField
forall b c a. (b -> c) -> (a -> b) -> a -> c
. VectorField
vf
          vy :: ScalarField
vy = Vec -> R
yComp (Vec -> R) -> VectorField -> ScalarField
forall b c a. (b -> c) -> (a -> b) -> a -> c
. VectorField
vf
          vz :: ScalarField
vz = Vec -> R
zComp (Vec -> R) -> VectorField -> ScalarField
forall b c a. (b -> c) -> (a -> b) -> a -> c
. VectorField
vf
          derivX :: ScalarField -> ScalarField
derivX = Vec -> ScalarField -> ScalarField
directionalDerivative (R
a R -> Vec -> Vec
*^ Vec
iHat)
          derivY :: ScalarField -> ScalarField
derivY = Vec -> ScalarField -> ScalarField
directionalDerivative (R
a R -> Vec -> Vec
*^ Vec
jHat)
          derivZ :: ScalarField -> ScalarField
derivZ = Vec -> ScalarField -> ScalarField
directionalDerivative (R
a R -> Vec -> Vec
*^ Vec
kHat)
      in      (ScalarField -> ScalarField
derivY ScalarField
vz Position
r R -> R -> R
forall a. Num a => a -> a -> a
- ScalarField -> ScalarField
derivZ ScalarField
vy Position
r) R -> Vec -> Vec
*^ Vec
iHat
          Vec -> Vec -> Vec
^+^ (ScalarField -> ScalarField
derivZ ScalarField
vx Position
r R -> R -> R
forall a. Num a => a -> a -> a
- ScalarField -> ScalarField
derivX ScalarField
vz Position
r) R -> Vec -> Vec
*^ Vec
jHat
          Vec -> Vec -> Vec
^+^ (ScalarField -> ScalarField
derivX ScalarField
vy Position
r R -> R -> R
forall a. Num a => a -> a -> a
- ScalarField -> ScalarField
derivY ScalarField
vx Position
r) R -> Vec -> Vec
*^ Vec
kHat

type FieldState = (R            -- time t
                  ,VectorField  -- electric field E
                  ,VectorField  -- magnetic field B
                  )

maxwellUpdate :: R                   -- dx
              -> R                   -- dt
              -> (R -> VectorField)  -- J
              -> FieldState -> FieldState
maxwellUpdate :: R -> R -> (R -> VectorField) -> FieldState -> FieldState
maxwellUpdate R
dx R
dt R -> VectorField
j (R
t,VectorField
eF,VectorField
bF)
    = let t' :: R
t'    = R
t R -> R -> R
forall a. Num a => a -> a -> a
+ R
dt
          eF' :: VectorField
eF' Position
r = VectorField
eF Position
r Vec -> Vec -> Vec
^+^ R
cSIR -> R -> R
forall a. Floating a => a -> a -> a
**R
2 R -> Vec -> Vec
*^ R
dt R -> Vec -> Vec
*^ (R -> VectorField -> VectorField
curl R
dx VectorField
bF Position
r Vec -> Vec -> Vec
^-^ R
mu0 R -> Vec -> Vec
*^ R -> VectorField
j R
t Position
r)
          bF' :: VectorField
bF' Position
r = VectorField
bF Position
r Vec -> Vec -> Vec
^-^           R
dt R -> Vec -> Vec
*^  R -> VectorField -> VectorField
curl R
dx VectorField
eF Position
r
      in (R
t',VectorField
eF',VectorField
bF')

maxwellEvolve :: R                   -- dx
              -> R                   -- dt
              -> (R -> VectorField)  -- J
              -> FieldState -> [FieldState]
maxwellEvolve :: R -> R -> (R -> VectorField) -> FieldState -> [FieldState]
maxwellEvolve R
dx R
dt R -> VectorField
j FieldState
st0 = (FieldState -> FieldState) -> FieldState -> [FieldState]
forall a. (a -> a) -> a -> [a]
iterate (R -> R -> (R -> VectorField) -> FieldState -> FieldState
maxwellUpdate R
dx R
dt R -> VectorField
j) FieldState
st0

exLocs, eyLocs, ezLocs, bxLocs, byLocs, bzLocs :: [(Int,Int,Int)]
exLocs :: [(Int, Int, Int)]
exLocs = [(Int
nx,Int
ny,Int
nz) | Int
nx <- [Int]
odds , Int
ny <- [Int]
evens, Int
nz <- [Int]
evens]
eyLocs :: [(Int, Int, Int)]
eyLocs = [(Int
nx,Int
ny,Int
nz) | Int
nx <- [Int]
evens, Int
ny <- [Int]
odds , Int
nz <- [Int]
evens]
ezLocs :: [(Int, Int, Int)]
ezLocs = [(Int
nx,Int
ny,Int
nz) | Int
nx <- [Int]
evens, Int
ny <- [Int]
evens, Int
nz <- [Int]
odds ]
bxLocs :: [(Int, Int, Int)]
bxLocs = [(Int
nx,Int
ny,Int
nz) | Int
nx <- [Int]
evens, Int
ny <- [Int]
odds , Int
nz <- [Int]
odds ]
byLocs :: [(Int, Int, Int)]
byLocs = [(Int
nx,Int
ny,Int
nz) | Int
nx <- [Int]
odds , Int
ny <- [Int]
evens, Int
nz <- [Int]
odds ]
bzLocs :: [(Int, Int, Int)]
bzLocs = [(Int
nx,Int
ny,Int
nz) | Int
nx <- [Int]
odds , Int
ny <- [Int]
odds , Int
nz <- [Int]
evens]

spaceStepsCE :: Int
spaceStepsCE :: Int
spaceStepsCE = Int
40

hiEven :: Int
hiEven :: Int
hiEven =  Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
spaceStepsCE

evens :: [Int]
evens :: [Int]
evens = [-Int
hiEven, -Int
hiEven Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2 .. Int
hiEven]

odds :: [Int]
odds :: [Int]
odds = [-Int
hiEven Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1, -Int
hiEven Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
3 .. Int
hiEven Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]

data StateFDTD = StateFDTD {StateFDTD -> R
timeFDTD :: R
                           ,StateFDTD -> R
stepX    :: R
                           ,StateFDTD -> R
stepY    :: R
                           ,StateFDTD -> R
stepZ    :: R
                           ,StateFDTD -> Map (Int, Int, Int) R
eField   :: M.Map (Int,Int,Int) R
                           ,StateFDTD -> Map (Int, Int, Int) R
bField   :: M.Map (Int,Int,Int) R
                           } deriving Int -> StateFDTD -> ShowS
[StateFDTD] -> ShowS
StateFDTD -> String
(Int -> StateFDTD -> ShowS)
-> (StateFDTD -> String)
-> ([StateFDTD] -> ShowS)
-> Show StateFDTD
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> StateFDTD -> ShowS
showsPrec :: Int -> StateFDTD -> ShowS
$cshow :: StateFDTD -> String
show :: StateFDTD -> String
$cshowList :: [StateFDTD] -> ShowS
showList :: [StateFDTD] -> ShowS
Show

initialStateFDTD :: R -> StateFDTD
initialStateFDTD :: R -> StateFDTD
initialStateFDTD R
spatialStep
    = StateFDTD {timeFDTD :: R
timeFDTD  = R
0
                ,stepX :: R
stepX = R
spatialStep
                ,stepY :: R
stepY = R
spatialStep
                ,stepZ :: R
stepZ = R
spatialStep
                ,eField :: Map (Int, Int, Int) R
eField = [((Int, Int, Int), R)] -> Map (Int, Int, Int) R
forall k a. Ord k => [(k, a)] -> Map k a
M.fromList [((Int, Int, Int)
loc,R
0) | (Int, Int, Int)
loc <- [(Int, Int, Int)]
exLocs[(Int, Int, Int)] -> [(Int, Int, Int)] -> [(Int, Int, Int)]
forall a. [a] -> [a] -> [a]
++[(Int, Int, Int)]
eyLocs[(Int, Int, Int)] -> [(Int, Int, Int)] -> [(Int, Int, Int)]
forall a. [a] -> [a] -> [a]
++[(Int, Int, Int)]
ezLocs]
                ,bField :: Map (Int, Int, Int) R
bField = [((Int, Int, Int), R)] -> Map (Int, Int, Int) R
forall k a. Ord k => [(k, a)] -> Map k a
M.fromList [((Int, Int, Int)
loc,R
0) | (Int, Int, Int)
loc <- [(Int, Int, Int)]
bxLocs[(Int, Int, Int)] -> [(Int, Int, Int)] -> [(Int, Int, Int)]
forall a. [a] -> [a] -> [a]
++[(Int, Int, Int)]
byLocs[(Int, Int, Int)] -> [(Int, Int, Int)] -> [(Int, Int, Int)]
forall a. [a] -> [a] -> [a]
++[(Int, Int, Int)]
bzLocs]
                }

lookupAZ :: Ord k => k -> M.Map k R -> R
lookupAZ :: forall k. Ord k => k -> Map k R -> R
lookupAZ k
key Map k R
m = case k -> Map k R -> Maybe R
forall k a. Ord k => k -> Map k a -> Maybe a
M.lookup k
key Map k R
m of
                     Maybe R
Nothing -> R
0
                     Just R
x  -> R
x

partialX,partialY,partialZ :: R -> M.Map (Int,Int,Int) R -> (Int,Int,Int) -> R
partialX :: R -> Map (Int, Int, Int) R -> (Int, Int, Int) -> R
partialX R
dx Map (Int, Int, Int) R
m (Int
i,Int
j,Int
k) = ((Int, Int, Int) -> Map (Int, Int, Int) R -> R
forall k. Ord k => k -> Map k R -> R
lookupAZ (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1,Int
j,Int
k) Map (Int, Int, Int) R
m R -> R -> R
forall a. Num a => a -> a -> a
- (Int, Int, Int) -> Map (Int, Int, Int) R -> R
forall k. Ord k => k -> Map k R -> R
lookupAZ (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1,Int
j,Int
k) Map (Int, Int, Int) R
m) R -> R -> R
forall a. Fractional a => a -> a -> a
/ R
dx
partialY :: R -> Map (Int, Int, Int) R -> (Int, Int, Int) -> R
partialY R
dy Map (Int, Int, Int) R
m (Int
i,Int
j,Int
k) = ((Int, Int, Int) -> Map (Int, Int, Int) R -> R
forall k. Ord k => k -> Map k R -> R
lookupAZ (Int
i,Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1,Int
k) Map (Int, Int, Int) R
m R -> R -> R
forall a. Num a => a -> a -> a
- (Int, Int, Int) -> Map (Int, Int, Int) R -> R
forall k. Ord k => k -> Map k R -> R
lookupAZ (Int
i,Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1,Int
k) Map (Int, Int, Int) R
m) R -> R -> R
forall a. Fractional a => a -> a -> a
/ R
dy
partialZ :: R -> Map (Int, Int, Int) R -> (Int, Int, Int) -> R
partialZ R
dz Map (Int, Int, Int) R
m (Int
i,Int
j,Int
k) = ((Int, Int, Int) -> Map (Int, Int, Int) R -> R
forall k. Ord k => k -> Map k R -> R
lookupAZ (Int
i,Int
j,Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Map (Int, Int, Int) R
m R -> R -> R
forall a. Num a => a -> a -> a
- (Int, Int, Int) -> Map (Int, Int, Int) R -> R
forall k. Ord k => k -> Map k R -> R
lookupAZ (Int
i,Int
j,Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) Map (Int, Int, Int) R
m) R -> R -> R
forall a. Fractional a => a -> a -> a
/ R
dz

curlEx,curlEy,curlEz,curlBx,curlBy,curlBz :: StateFDTD -> (Int,Int,Int) -> R
curlBx :: StateFDTD -> (Int, Int, Int) -> R
curlBx (StateFDTD R
_ R
_ R
dy R
dz Map (Int, Int, Int) R
_ Map (Int, Int, Int) R
b) (Int, Int, Int)
loc = R -> Map (Int, Int, Int) R -> (Int, Int, Int) -> R
partialY R
dy Map (Int, Int, Int) R
b (Int, Int, Int)
loc R -> R -> R
forall a. Num a => a -> a -> a
- R -> Map (Int, Int, Int) R -> (Int, Int, Int) -> R
partialZ R
dz Map (Int, Int, Int) R
b (Int, Int, Int)
loc
curlBy :: StateFDTD -> (Int, Int, Int) -> R
curlBy (StateFDTD R
_ R
dx R
_ R
dz Map (Int, Int, Int) R
_ Map (Int, Int, Int) R
b) (Int, Int, Int)
loc = R -> Map (Int, Int, Int) R -> (Int, Int, Int) -> R
partialZ R
dz Map (Int, Int, Int) R
b (Int, Int, Int)
loc R -> R -> R
forall a. Num a => a -> a -> a
- R -> Map (Int, Int, Int) R -> (Int, Int, Int) -> R
partialX R
dx Map (Int, Int, Int) R
b (Int, Int, Int)
loc
curlBz :: StateFDTD -> (Int, Int, Int) -> R
curlBz (StateFDTD R
_ R
dx R
dy R
_ Map (Int, Int, Int) R
_ Map (Int, Int, Int) R
b) (Int, Int, Int)
loc = R -> Map (Int, Int, Int) R -> (Int, Int, Int) -> R
partialX R
dx Map (Int, Int, Int) R
b (Int, Int, Int)
loc R -> R -> R
forall a. Num a => a -> a -> a
- R -> Map (Int, Int, Int) R -> (Int, Int, Int) -> R
partialY R
dy Map (Int, Int, Int) R
b (Int, Int, Int)
loc
curlEx :: StateFDTD -> (Int, Int, Int) -> R
curlEx (StateFDTD R
_ R
_ R
dy R
dz Map (Int, Int, Int) R
e Map (Int, Int, Int) R
_) (Int, Int, Int)
loc = R -> Map (Int, Int, Int) R -> (Int, Int, Int) -> R
partialY R
dy Map (Int, Int, Int) R
e (Int, Int, Int)
loc R -> R -> R
forall a. Num a => a -> a -> a
- R -> Map (Int, Int, Int) R -> (Int, Int, Int) -> R
partialZ R
dz Map (Int, Int, Int) R
e (Int, Int, Int)
loc
curlEy :: StateFDTD -> (Int, Int, Int) -> R
curlEy (StateFDTD R
_ R
dx R
_ R
dz Map (Int, Int, Int) R
e Map (Int, Int, Int) R
_) (Int, Int, Int)
loc = R -> Map (Int, Int, Int) R -> (Int, Int, Int) -> R
partialZ R
dz Map (Int, Int, Int) R
e (Int, Int, Int)
loc R -> R -> R
forall a. Num a => a -> a -> a
- R -> Map (Int, Int, Int) R -> (Int, Int, Int) -> R
partialX R
dx Map (Int, Int, Int) R
e (Int, Int, Int)
loc
curlEz :: StateFDTD -> (Int, Int, Int) -> R
curlEz (StateFDTD R
_ R
dx R
dy R
_ Map (Int, Int, Int) R
e Map (Int, Int, Int) R
_) (Int, Int, Int)
loc = R -> Map (Int, Int, Int) R -> (Int, Int, Int) -> R
partialX R
dx Map (Int, Int, Int) R
e (Int, Int, Int)
loc R -> R -> R
forall a. Num a => a -> a -> a
- R -> Map (Int, Int, Int) R -> (Int, Int, Int) -> R
partialY R
dy Map (Int, Int, Int) R
e (Int, Int, Int)
loc

stateUpdate :: R                   -- dt
            -> (R -> VectorField)  -- current density J
            -> StateFDTD -> StateFDTD
stateUpdate :: R -> (R -> VectorField) -> StateFDTD -> StateFDTD
stateUpdate R
dt R -> VectorField
j st0 :: StateFDTD
st0@(StateFDTD R
t R
_dx R
_dy R
_dz Map (Int, Int, Int) R
_e Map (Int, Int, Int) R
_b)
    = let st1 :: StateFDTD
st1 = R -> VectorField -> StateFDTD -> StateFDTD
updateE R
dt (R -> VectorField
j R
t) StateFDTD
st0
          st2 :: StateFDTD
st2 = R -> StateFDTD -> StateFDTD
updateB R
dt StateFDTD
st1
      in StateFDTD
st2

updateE :: R            -- time step dt
        -> VectorField  -- current density J
        -> StateFDTD -> StateFDTD
updateE :: R -> VectorField -> StateFDTD -> StateFDTD
updateE R
dt VectorField
jVF StateFDTD
st
    = StateFDTD
st { timeFDTD = timeFDTD st + dt / 2
         , eField   = M.mapWithKey (updateEOneLoc dt jVF st) (eField st) }

updateB :: R -> StateFDTD -> StateFDTD
updateB :: R -> StateFDTD -> StateFDTD
updateB R
dt StateFDTD
st
    = StateFDTD
st { timeFDTD = timeFDTD st + dt / 2
         , bField   = M.mapWithKey (updateBOneLoc dt st) (bField st) }

updateEOneLoc :: R -> VectorField -> StateFDTD -> (Int,Int,Int) -> R -> R
updateEOneLoc :: R -> VectorField -> StateFDTD -> (Int, Int, Int) -> R -> R
updateEOneLoc R
dt VectorField
jVF StateFDTD
st (Int
nx,Int
ny,Int
nz) R
ec
    = let r :: Position
r = R -> R -> R -> Position
cart (Int -> R
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
nx R -> R -> R
forall a. Num a => a -> a -> a
* StateFDTD -> R
stepX StateFDTD
st R -> R -> R
forall a. Fractional a => a -> a -> a
/ R
2)
                   (Int -> R
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
ny R -> R -> R
forall a. Num a => a -> a -> a
* StateFDTD -> R
stepY StateFDTD
st R -> R -> R
forall a. Fractional a => a -> a -> a
/ R
2)
                   (Int -> R
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
nz R -> R -> R
forall a. Num a => a -> a -> a
* StateFDTD -> R
stepZ StateFDTD
st R -> R -> R
forall a. Fractional a => a -> a -> a
/ R
2)
          Vec R
jx R
jy R
jz = VectorField
jVF Position
r
      in case (Int -> Bool
forall a. Integral a => a -> Bool
odd Int
nx, Int -> Bool
forall a. Integral a => a -> Bool
odd Int
ny, Int -> Bool
forall a. Integral a => a -> Bool
odd Int
nz) of
           (Bool
True , Bool
False, Bool
False)
               -> R
ec R -> R -> R
forall a. Num a => a -> a -> a
+ R
cSIR -> R -> R
forall a. Floating a => a -> a -> a
**R
2 R -> R -> R
forall a. Num a => a -> a -> a
* (StateFDTD -> (Int, Int, Int) -> R
curlBx StateFDTD
st (Int
nx,Int
ny,Int
nz) R -> R -> R
forall a. Num a => a -> a -> a
- R
mu0 R -> R -> R
forall a. Num a => a -> a -> a
* R
jx) R -> R -> R
forall a. Num a => a -> a -> a
* R
dt  -- Ex
           (Bool
False, Bool
True , Bool
False)
               -> R
ec R -> R -> R
forall a. Num a => a -> a -> a
+ R
cSIR -> R -> R
forall a. Floating a => a -> a -> a
**R
2 R -> R -> R
forall a. Num a => a -> a -> a
* (StateFDTD -> (Int, Int, Int) -> R
curlBy StateFDTD
st (Int
nx,Int
ny,Int
nz) R -> R -> R
forall a. Num a => a -> a -> a
- R
mu0 R -> R -> R
forall a. Num a => a -> a -> a
* R
jy) R -> R -> R
forall a. Num a => a -> a -> a
* R
dt  -- Ey
           (Bool
False, Bool
False, Bool
True )
               -> R
ec R -> R -> R
forall a. Num a => a -> a -> a
+ R
cSIR -> R -> R
forall a. Floating a => a -> a -> a
**R
2 R -> R -> R
forall a. Num a => a -> a -> a
* (StateFDTD -> (Int, Int, Int) -> R
curlBz StateFDTD
st (Int
nx,Int
ny,Int
nz) R -> R -> R
forall a. Num a => a -> a -> a
- R
mu0 R -> R -> R
forall a. Num a => a -> a -> a
* R
jz) R -> R -> R
forall a. Num a => a -> a -> a
* R
dt  -- Ez
           (Bool, Bool, Bool)
_ -> String -> R
forall a. HasCallStack => String -> a
error String
"updateEOneLoc passed bad indices"

updateBOneLoc :: R -> StateFDTD -> (Int,Int,Int) -> R -> R
updateBOneLoc :: R -> StateFDTD -> (Int, Int, Int) -> R -> R
updateBOneLoc R
dt StateFDTD
st (Int
nx,Int
ny,Int
nz) R
bc
    = case (Int -> Bool
forall a. Integral a => a -> Bool
odd Int
nx, Int -> Bool
forall a. Integral a => a -> Bool
odd Int
ny, Int -> Bool
forall a. Integral a => a -> Bool
odd Int
nz) of
        (Bool
False, Bool
True , Bool
True ) -> R
bc R -> R -> R
forall a. Num a => a -> a -> a
- StateFDTD -> (Int, Int, Int) -> R
curlEx StateFDTD
st (Int
nx,Int
ny,Int
nz) R -> R -> R
forall a. Num a => a -> a -> a
* R
dt  -- Bx
        (Bool
True , Bool
False, Bool
True ) -> R
bc R -> R -> R
forall a. Num a => a -> a -> a
- StateFDTD -> (Int, Int, Int) -> R
curlEy StateFDTD
st (Int
nx,Int
ny,Int
nz) R -> R -> R
forall a. Num a => a -> a -> a
* R
dt  -- By
        (Bool
True , Bool
True , Bool
False) -> R
bc R -> R -> R
forall a. Num a => a -> a -> a
- StateFDTD -> (Int, Int, Int) -> R
curlEz StateFDTD
st (Int
nx,Int
ny,Int
nz) R -> R -> R
forall a. Num a => a -> a -> a
* R
dt  -- Bz
        (Bool, Bool, Bool)
_ -> String -> R
forall a. HasCallStack => String -> a
error String
"updateBOneLoc passed bad indices"

jGaussian :: R -> VectorField
jGaussian :: R -> VectorField
jGaussian R
t Position
r
    = let wavelength :: R
wavelength = R
1.08             -- meters
          frequency :: R
frequency = R
cSI R -> R -> R
forall a. Fractional a => a -> a -> a
/ R
wavelength  -- Hz
          j0 :: R
j0 = R
77.5                     -- A/m^2
          l :: R
l = R
0.108                     -- meters
          rMag :: R
rMag = Vec -> R
magnitude (VectorField
rVF Position
r)      -- meters
      in R
j0 R -> Vec -> Vec
*^ R -> R
forall a. Floating a => a -> a
exp (-R
rMagR -> R -> R
forall a. Floating a => a -> a -> a
**R
2 R -> R -> R
forall a. Fractional a => a -> a -> a
/ R
lR -> R -> R
forall a. Floating a => a -> a -> a
**R
2) R -> Vec -> Vec
*^ R -> R
forall a. Floating a => a -> a
cos (R
2R -> R -> R
forall a. Num a => a -> a -> a
*R
forall a. Floating a => a
piR -> R -> R
forall a. Num a => a -> a -> a
*R
frequencyR -> R -> R
forall a. Num a => a -> a -> a
*R
t) R -> Vec -> Vec
*^ Vec
kHat

getAverage :: (Int,Int,Int)  -- (even,even,even) or (odd,odd,odd)
           -> M.Map (Int,Int,Int) R
           -> Vec
getAverage :: (Int, Int, Int) -> Map (Int, Int, Int) R -> Vec
getAverage (Int
i,Int
j,Int
k) Map (Int, Int, Int) R
m
    = let vXl :: R
vXl = (Int, Int, Int) -> Map (Int, Int, Int) R -> R
forall k. Ord k => k -> Map k R -> R
lookupAZ (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1,Int
j  ,Int
k  ) Map (Int, Int, Int) R
m
          vYl :: R
vYl = (Int, Int, Int) -> Map (Int, Int, Int) R -> R
forall k. Ord k => k -> Map k R -> R
lookupAZ (Int
i  ,Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1,Int
k  ) Map (Int, Int, Int) R
m
          vZl :: R
vZl = (Int, Int, Int) -> Map (Int, Int, Int) R -> R
forall k. Ord k => k -> Map k R -> R
lookupAZ (Int
i  ,Int
j  ,Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) Map (Int, Int, Int) R
m
          vXr :: R
vXr = (Int, Int, Int) -> Map (Int, Int, Int) R -> R
forall k. Ord k => k -> Map k R -> R
lookupAZ (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1,Int
j  ,Int
k  ) Map (Int, Int, Int) R
m
          vYr :: R
vYr = (Int, Int, Int) -> Map (Int, Int, Int) R -> R
forall k. Ord k => k -> Map k R -> R
lookupAZ (Int
i  ,Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1,Int
k  ) Map (Int, Int, Int) R
m
          vZr :: R
vZr = (Int, Int, Int) -> Map (Int, Int, Int) R -> R
forall k. Ord k => k -> Map k R -> R
lookupAZ (Int
i  ,Int
j  ,Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Map (Int, Int, Int) R
m
      in R -> R -> R -> Vec
vec ((R
vXlR -> R -> R
forall a. Num a => a -> a -> a
+R
vXr)R -> R -> R
forall a. Fractional a => a -> a -> a
/R
2) ((R
vYlR -> R -> R
forall a. Num a => a -> a -> a
+R
vYr)R -> R -> R
forall a. Fractional a => a -> a -> a
/R
2) ((R
vZlR -> R -> R
forall a. Num a => a -> a -> a
+R
vZr)R -> R -> R
forall a. Fractional a => a -> a -> a
/R
2)