{-# OPTIONS -Wall #-}

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

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

module LPFPCore.ElectricField where

import LPFPCore.SimpleVec
    ( R, Vec, (^+^), (^-^), (*^), (^*), (^/), (<.>), (><)
    , sumV, magnitude, vec )
import LPFPCore.CoordinateSystems
    ( Position, ScalarField, VectorField
    , displacement, shiftPosition, addVectorFields
    , origin, rVF )
import LPFPCore.Geometry ( Curve(..), Surface(..), Volume(..) )
import LPFPCore.Charge
    ( Charge, ChargeDistribution(..)
    , diskCap, simpleDipole, lineDipole )

epsilon0 :: R
epsilon0 :: R
epsilon0 = R
1forall a. Fractional a => a -> a -> a
/(R
mu0 forall a. Num a => a -> a -> a
* R
cSIforall a. Floating a => a -> a -> a
**R
2)

cSI :: R
cSI :: R
cSI = R
299792458  -- m/s

mu0 :: R
mu0 :: R
mu0 = R
4e-7 forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi  -- N/A^2

eFieldFromPointCharge
    :: Charge       -- in Coulombs
    -> Position     -- of point charge (in m)
    -> VectorField  -- electric field (in V/m)
eFieldFromPointCharge :: R -> Position -> VectorField
eFieldFromPointCharge R
q1 Position
r1 Position
r
    = let k :: R
k = R
1 forall a. Fractional a => a -> a -> a
/ (R
4 forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* R
epsilon0)
          d :: Displacement
d = Position -> VectorField
displacement Position
r1 Position
r
      in (R
k forall a. Num a => a -> a -> a
* R
q1) R -> Displacement -> Displacement
*^ Displacement
d Displacement -> R -> Displacement
^/ Displacement -> R
magnitude Displacement
d forall a. Floating a => a -> a -> a
** R
3

eField :: ChargeDistribution -> VectorField
eField :: ChargeDistribution -> VectorField
eField (PointCharge   R
q   Position
r) = R -> Position -> VectorField
eFieldFromPointCharge   R
q   Position
r
eField (LineCharge    ScalarField
lam Curve
c) = ScalarField -> Curve -> VectorField
eFieldFromLineCharge    ScalarField
lam Curve
c
eField (SurfaceCharge ScalarField
sig Surface
s) = ScalarField -> Surface -> VectorField
eFieldFromSurfaceCharge ScalarField
sig Surface
s
eField (VolumeCharge  ScalarField
rho Volume
v) = ScalarField -> Volume -> VectorField
eFieldFromVolumeCharge  ScalarField
rho Volume
v
eField (MultipleCharges [ChargeDistribution]
cds) = [VectorField] -> VectorField
addVectorFields forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map ChargeDistribution -> VectorField
eField [ChargeDistribution]
cds

simpleDipoleSodiumChloride :: ChargeDistribution
simpleDipoleSodiumChloride :: ChargeDistribution
simpleDipoleSodiumChloride = Displacement -> R -> ChargeDistribution
simpleDipole (R -> R -> R -> Displacement
vec R
0 R
0 R
2.99e-29) R
2.36e-10

eFieldSodiumChloride :: VectorField
eFieldSodiumChloride :: VectorField
eFieldSodiumChloride = ChargeDistribution -> VectorField
eField ChargeDistribution
simpleDipoleSodiumChloride

eFieldIdealDipole :: Vec          -- electric dipole moment
                  -> VectorField  -- electric field
eFieldIdealDipole :: Displacement -> VectorField
eFieldIdealDipole Displacement
p Position
r
    = let k :: R
k = R
1 forall a. Fractional a => a -> a -> a
/ (R
4 forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* R
epsilon0)  -- SI units
          rMag :: R
rMag = Displacement -> R
magnitude (VectorField
rVF Position
r)
          rUnit :: Displacement
rUnit = VectorField
rVF Position
r Displacement -> R -> Displacement
^/ R
rMag
      in R
k R -> Displacement -> Displacement
*^ (R
1 forall a. Fractional a => a -> a -> a
/ R
rMagforall a. Floating a => a -> a -> a
**R
3) R -> Displacement -> Displacement
*^ (R
3 R -> Displacement -> Displacement
*^ (Displacement
p Displacement -> Displacement -> R
<.> Displacement
rUnit) R -> Displacement -> Displacement
*^ Displacement
rUnit Displacement -> Displacement -> Displacement
^-^ Displacement
p)

type VectorLineIntegral = VectorField -> Curve -> Vec

type CurveApprox = Curve -> [(Position,Vec)]

vectorLineIntegral :: CurveApprox -> VectorField -> Curve -> Vec
vectorLineIntegral :: CurveApprox -> VectorField -> Curve -> Displacement
vectorLineIntegral CurveApprox
approx VectorField
vF Curve
c
    = [Displacement] -> Displacement
sumV [VectorField
vF Position
r' Displacement -> R -> Displacement
^* Displacement -> R
magnitude Displacement
dl' | (Position
r',Displacement
dl') <- CurveApprox
approx Curve
c]

eFieldFromLineCharge
    :: ScalarField  -- linear charge density lambda
    -> Curve        -- geometry of the line charge
    -> VectorField  -- electric field (in V/m)
eFieldFromLineCharge :: ScalarField -> Curve -> VectorField
eFieldFromLineCharge ScalarField
lambda Curve
c Position
r
    = let k :: R
k = R
1 forall a. Fractional a => a -> a -> a
/ (R
4 forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* R
epsilon0)
          integrand :: VectorField
integrand Position
r' = ScalarField
lambda Position
r' R -> Displacement -> Displacement
*^ Displacement
d Displacement -> R -> Displacement
^/ Displacement -> R
magnitude Displacement
d forall a. Floating a => a -> a -> a
** R
3
              where d :: Displacement
d = Position -> VectorField
displacement Position
r' Position
r
      in R
k R -> Displacement -> Displacement
*^ CurveApprox -> VectorField -> Curve -> Displacement
vectorLineIntegral (Int -> CurveApprox
curveSample Int
1000) VectorField
integrand Curve
c

lineDipoleSodiumChloride :: ChargeDistribution
lineDipoleSodiumChloride :: ChargeDistribution
lineDipoleSodiumChloride = Displacement -> R -> ChargeDistribution
lineDipole (R -> R -> R -> Displacement
vec R
0 R
0 R
2.99e-29) R
2.36e-10

eFieldLineDipole :: VectorField
eFieldLineDipole :: VectorField
eFieldLineDipole = ChargeDistribution -> VectorField
eField ChargeDistribution
lineDipoleSodiumChloride

type VectorSurfaceIntegral = VectorField -> Surface -> Vec

type SurfaceApprox = Surface -> [(Position,Vec)]

vectorSurfaceIntegral :: SurfaceApprox -> VectorField -> Surface -> Vec
vectorSurfaceIntegral :: SurfaceApprox -> VectorField -> Surface -> Displacement
vectorSurfaceIntegral SurfaceApprox
approx VectorField
vF Surface
s
    = [Displacement] -> Displacement
sumV [VectorField
vF Position
r' Displacement -> R -> Displacement
^* Displacement -> R
magnitude Displacement
da' | (Position
r',Displacement
da') <- SurfaceApprox
approx Surface
s]

eFieldFromSurfaceCharge
    :: ScalarField  -- surface charge density sigma
    -> Surface      -- geometry of the surface charge
    -> VectorField  -- electric field (in V/m)
eFieldFromSurfaceCharge :: ScalarField -> Surface -> VectorField
eFieldFromSurfaceCharge ScalarField
sigma Surface
s Position
r
    = let k :: R
k = R
1 forall a. Fractional a => a -> a -> a
/ (R
4 forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* R
epsilon0)
          integrand :: VectorField
integrand Position
r' = ScalarField
sigma Position
r' R -> Displacement -> Displacement
*^ Displacement
d Displacement -> R -> Displacement
^/ Displacement -> R
magnitude Displacement
d forall a. Floating a => a -> a -> a
** R
3
              where d :: Displacement
d = Position -> VectorField
displacement Position
r' Position
r
      in R
k R -> Displacement -> Displacement
*^ SurfaceApprox -> VectorField -> Surface -> Displacement
vectorSurfaceIntegral (Int -> SurfaceApprox
surfaceSample Int
200) VectorField
integrand Surface
s

eFieldDiskCap :: VectorField
eFieldDiskCap :: VectorField
eFieldDiskCap = ChargeDistribution -> VectorField
eField forall a b. (a -> b) -> a -> b
$ R -> R -> R -> ChargeDistribution
diskCap R
0.05 R
0.04 R
2e-8

type VectorVolumeIntegral = VectorField -> Volume -> Vec

type VolumeApprox = Volume -> [(Position,R)]

vectorVolumeIntegral :: VolumeApprox -> VectorField -> Volume -> Vec
vectorVolumeIntegral :: VolumeApprox -> VectorField -> Volume -> Displacement
vectorVolumeIntegral VolumeApprox
approx VectorField
vF Volume
vol
    = [Displacement] -> Displacement
sumV [VectorField
vF Position
r' Displacement -> R -> Displacement
^* R
dv' | (Position
r',R
dv') <- VolumeApprox
approx Volume
vol]

eFieldFromVolumeCharge
    :: ScalarField  -- volume charge density rho
    -> Volume       -- geometry of the volume charge
    -> VectorField  -- electric field (in V/m)
eFieldFromVolumeCharge :: ScalarField -> Volume -> VectorField
eFieldFromVolumeCharge ScalarField
rho Volume
v Position
r
    = let k :: R
k = R
1 forall a. Fractional a => a -> a -> a
/ (R
4 forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* R
epsilon0)
          integrand :: VectorField
integrand Position
r' = ScalarField
rho Position
r' R -> Displacement -> Displacement
*^ Displacement
d Displacement -> R -> Displacement
^/ Displacement -> R
magnitude Displacement
d forall a. Floating a => a -> a -> a
** R
3
              where d :: Displacement
d = Position -> VectorField
displacement Position
r' Position
r
      in R
k R -> Displacement -> Displacement
*^ VolumeApprox -> VectorField -> Volume -> Displacement
vectorVolumeIntegral (Int -> VolumeApprox
volumeSample Int
50) VectorField
integrand Volume
v

type ScalarLineIntegral = ScalarField -> Curve -> R

scalarLineIntegral :: CurveApprox -> ScalarField -> Curve -> R
scalarLineIntegral :: CurveApprox -> ScalarField -> Curve -> R
scalarLineIntegral CurveApprox
approx ScalarField
f Curve
c
    = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [ScalarField
f Position
r' forall a. Num a => a -> a -> a
* Displacement -> R
magnitude Displacement
dl' | (Position
r',Displacement
dl') <- CurveApprox
approx Curve
c]

type ScalarSurfaceIntegral = ScalarField -> Surface -> R

scalarSurfaceIntegral :: SurfaceApprox -> ScalarField -> Surface -> R
scalarSurfaceIntegral :: SurfaceApprox -> ScalarField -> Surface -> R
scalarSurfaceIntegral SurfaceApprox
approx ScalarField
f Surface
s
    = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [ScalarField
f Position
r' forall a. Num a => a -> a -> a
* Displacement -> R
magnitude Displacement
da' | (Position
r',Displacement
da') <- SurfaceApprox
approx Surface
s]

type ScalarVolumeIntegral = ScalarField -> Volume -> R

scalarVolumeIntegral :: VolumeApprox -> ScalarField -> Volume -> R
scalarVolumeIntegral :: VolumeApprox -> ScalarField -> Volume -> R
scalarVolumeIntegral VolumeApprox
approx ScalarField
f Volume
vol
    = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [ScalarField
f Position
r' forall a. Num a => a -> a -> a
* R
dv' | (Position
r',R
dv') <- VolumeApprox
approx Volume
vol]

curveSample :: Int -> Curve -> [(Position,Vec)]
curveSample :: Int -> CurveApprox
curveSample Int
n Curve
c
    = let segCent :: Segment -> Position
          segCent :: Segment -> Position
segCent (Position
p1,Position
p2) = Displacement -> Position -> Position
shiftPosition ((VectorField
rVF Position
p1 Displacement -> Displacement -> Displacement
^+^ VectorField
rVF Position
p2) Displacement -> R -> Displacement
^/ R
2) Position
origin
          segDisp :: Segment -> Vec
          segDisp :: Segment -> Displacement
segDisp = forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry Position -> VectorField
displacement
      in [(Segment -> Position
segCent Segment
seg, Segment -> Displacement
segDisp Segment
seg) | Segment
seg <- Int -> Curve -> [Segment]
segments Int
n Curve
c]

type Segment = (Position,Position)

segments :: Int -> Curve -> [Segment]
segments :: Int -> Curve -> [Segment]
segments Int
n (Curve R -> Position
g R
a R
b)
    = let ps :: [Position]
ps = forall a b. (a -> b) -> [a] -> [b]
map R -> Position
g forall a b. (a -> b) -> a -> b
$ Int -> R -> R -> [R]
linSpaced Int
n R
a R
b
      in forall a b. [a] -> [b] -> [(a, b)]
zip [Position]
ps (forall a. [a] -> [a]
tail [Position]
ps)

linSpaced :: Int -> R -> R -> [R]
linSpaced :: Int -> R -> R -> [R]
linSpaced Int
n R
x0 R
x1 = forall a. Int -> [a] -> [a]
take (Int
nforall a. Num a => a -> a -> a
+Int
1) [R
x0, R
x0forall a. Num a => a -> a -> a
+R
dx .. R
x1]
    where dx :: R
dx = (R
x1 forall a. Num a => a -> a -> a
- R
x0) forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n

surfaceSample :: Int -> Surface -> [(Position,Vec)]
surfaceSample :: Int -> SurfaceApprox
surfaceSample Int
n Surface
s = [(Triangle -> Position
triCenter Triangle
tri, Triangle -> Displacement
triArea Triangle
tri) | Triangle
tri <- Int -> Surface -> [Triangle]
triangles Int
n Surface
s]

data Triangle = Tri Position Position Position

triCenter :: Triangle -> Position
triCenter :: Triangle -> Position
triCenter (Tri Position
p1 Position
p2 Position
p3)
    = Displacement -> Position -> Position
shiftPosition ((VectorField
rVF Position
p1 Displacement -> Displacement -> Displacement
^+^ VectorField
rVF Position
p2 Displacement -> Displacement -> Displacement
^+^ VectorField
rVF Position
p3) Displacement -> R -> Displacement
^/ R
3) Position
origin

triArea :: Triangle -> Vec  -- vector area
triArea :: Triangle -> Displacement
triArea (Tri Position
p1 Position
p2 Position
p3) = R
0.5 R -> Displacement -> Displacement
*^ (Position -> VectorField
displacement Position
p1 Position
p2 Displacement -> Displacement -> Displacement
>< Position -> VectorField
displacement Position
p2 Position
p3)

triangles :: Int -> Surface -> [Triangle]
triangles :: Int -> Surface -> [Triangle]
triangles Int
n (Surface (R, R) -> Position
g R
sl R
su R -> R
tl R -> R
tu)
    = let sts :: [[(R, R)]]
sts = [[(R
s,R
t) | R
t <- Int -> R -> R -> [R]
linSpaced Int
n (R -> R
tl R
s) (R -> R
tu R
s)]
                     | R
s <- Int -> R -> R -> [R]
linSpaced Int
n R
sl R
su]
          stSquares :: [((R, R), (R, R), (R, R), (R, R))]
stSquares = [( [[(R, R)]]
sts forall a. [a] -> Int -> a
!! Int
j     forall a. [a] -> Int -> a
!! Int
k
                       , [[(R, R)]]
sts forall a. [a] -> Int -> a
!! (Int
jforall a. Num a => a -> a -> a
+Int
1) forall a. [a] -> Int -> a
!! Int
k
                       , [[(R, R)]]
sts forall a. [a] -> Int -> a
!! (Int
jforall a. Num a => a -> a -> a
+Int
1) forall a. [a] -> Int -> a
!! (Int
kforall a. Num a => a -> a -> a
+Int
1)
                       , [[(R, R)]]
sts forall a. [a] -> Int -> a
!! Int
j     forall a. [a] -> Int -> a
!! (Int
kforall a. Num a => a -> a -> a
+Int
1))
                      | Int
j <- [Int
0..Int
nforall a. Num a => a -> a -> a
-Int
1], Int
k <- [Int
0..Int
nforall a. Num a => a -> a -> a
-Int
1]]
          twoTriangles :: ((R, R), (R, R), (R, R), (R, R)) -> [Triangle]
twoTriangles ((R, R)
pp1,(R, R)
pp2,(R, R)
pp3,(R, R)
pp4)
              = [Position -> Position -> Position -> Triangle
Tri ((R, R) -> Position
g (R, R)
pp1) ((R, R) -> Position
g (R, R)
pp2) ((R, R) -> Position
g (R, R)
pp3),Position -> Position -> Position -> Triangle
Tri ((R, R) -> Position
g (R, R)
pp1) ((R, R) -> Position
g (R, R)
pp3) ((R, R) -> Position
g (R, R)
pp4)]
      in forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap ((R, R), (R, R), (R, R), (R, R)) -> [Triangle]
twoTriangles [((R, R), (R, R), (R, R), (R, R))]
stSquares

volumeSample :: Int -> Volume -> [(Position,R)]
volumeSample :: Int -> VolumeApprox
volumeSample Int
n Volume
v = [(Tet -> Position
tetCenter Tet
tet, Tet -> R
tetVolume Tet
tet) | Tet
tet <- Int -> Volume -> [Tet]
tetrahedrons Int
n Volume
v]

data Tet = Tet Position Position Position Position

tetCenter :: Tet -> Position
tetCenter :: Tet -> Position
tetCenter (Tet Position
p1 Position
p2 Position
p3 Position
p4)
    = Displacement -> Position -> Position
shiftPosition ((VectorField
rVF Position
p1 Displacement -> Displacement -> Displacement
^+^ VectorField
rVF Position
p2 Displacement -> Displacement -> Displacement
^+^ VectorField
rVF Position
p3 Displacement -> Displacement -> Displacement
^+^ VectorField
rVF Position
p4) Displacement -> R -> Displacement
^/ R
4) Position
origin

tetVolume :: Tet -> R
tetVolume :: Tet -> R
tetVolume (Tet Position
p1 Position
p2 Position
p3 Position
p4)
    = forall a. Num a => a -> a
abs forall a b. (a -> b) -> a -> b
$ (Displacement
d1 Displacement -> Displacement -> R
<.> (Displacement
d2 Displacement -> Displacement -> Displacement
>< Displacement
d3)) forall a. Fractional a => a -> a -> a
/ R
6
      where
        d1 :: Displacement
d1 = Position -> VectorField
displacement Position
p1 Position
p4
        d2 :: Displacement
d2 = Position -> VectorField
displacement Position
p2 Position
p4
        d3 :: Displacement
d3 = Position -> VectorField
displacement Position
p3 Position
p4

data ParamCube
    = PC { ParamCube -> (R, R, R)
v000 :: (R,R,R)
         , ParamCube -> (R, R, R)
v001 :: (R,R,R)
         , ParamCube -> (R, R, R)
v010 :: (R,R,R)
         , ParamCube -> (R, R, R)
v011 :: (R,R,R)
         , ParamCube -> (R, R, R)
v100 :: (R,R,R)
         , ParamCube -> (R, R, R)
v101 :: (R,R,R)
         , ParamCube -> (R, R, R)
v110 :: (R,R,R)
         , ParamCube -> (R, R, R)
v111 :: (R,R,R)
         }

tetrahedrons :: Int -> Volume -> [Tet]
tetrahedrons :: Int -> Volume -> [Tet]
tetrahedrons Int
n (Volume (R, R, R) -> Position
g R
sl R
su R -> R
tl R -> R
tu R -> R -> R
ul R -> R -> R
uu)
    = let stus :: [[[(R, R, R)]]]
stus = [[[(R
s,R
t,R
u) | R
u <- Int -> R -> R -> [R]
linSpaced Int
n (R -> R -> R
ul R
s R
t) (R -> R -> R
uu R
s R
t)]
                            | R
t <- Int -> R -> R -> [R]
linSpaced Int
n (R -> R
tl R
s) (R -> R
tu R
s)]
                            | R
s <- Int -> R -> R -> [R]
linSpaced Int
n R
sl R
su]
          stCubes :: [ParamCube]
stCubes = [(R, R, R)
-> (R, R, R)
-> (R, R, R)
-> (R, R, R)
-> (R, R, R)
-> (R, R, R)
-> (R, R, R)
-> (R, R, R)
-> ParamCube
PC ([[[(R, R, R)]]]
stus forall a. [a] -> Int -> a
!!  Int
j    forall a. [a] -> Int -> a
!!  Int
k    forall a. [a] -> Int -> a
!!  Int
l   )
                        ([[[(R, R, R)]]]
stus forall a. [a] -> Int -> a
!!  Int
j    forall a. [a] -> Int -> a
!!  Int
k    forall a. [a] -> Int -> a
!! (Int
lforall a. Num a => a -> a -> a
+Int
1))
                        ([[[(R, R, R)]]]
stus forall a. [a] -> Int -> a
!!  Int
j    forall a. [a] -> Int -> a
!! (Int
kforall a. Num a => a -> a -> a
+Int
1) forall a. [a] -> Int -> a
!!  Int
l   )
                        ([[[(R, R, R)]]]
stus forall a. [a] -> Int -> a
!!  Int
j    forall a. [a] -> Int -> a
!! (Int
kforall a. Num a => a -> a -> a
+Int
1) forall a. [a] -> Int -> a
!! (Int
lforall a. Num a => a -> a -> a
+Int
1))
                        ([[[(R, R, R)]]]
stus forall a. [a] -> Int -> a
!! (Int
jforall a. Num a => a -> a -> a
+Int
1) forall a. [a] -> Int -> a
!!  Int
k    forall a. [a] -> Int -> a
!!  Int
l   )
                        ([[[(R, R, R)]]]
stus forall a. [a] -> Int -> a
!! (Int
jforall a. Num a => a -> a -> a
+Int
1) forall a. [a] -> Int -> a
!!  Int
k    forall a. [a] -> Int -> a
!! (Int
lforall a. Num a => a -> a -> a
+Int
1))
                        ([[[(R, R, R)]]]
stus forall a. [a] -> Int -> a
!! (Int
jforall a. Num a => a -> a -> a
+Int
1) forall a. [a] -> Int -> a
!! (Int
kforall a. Num a => a -> a -> a
+Int
1) forall a. [a] -> Int -> a
!!  Int
l   )
                        ([[[(R, R, R)]]]
stus forall a. [a] -> Int -> a
!! (Int
jforall a. Num a => a -> a -> a
+Int
1) forall a. [a] -> Int -> a
!! (Int
kforall a. Num a => a -> a -> a
+Int
1) forall a. [a] -> Int -> a
!! (Int
lforall a. Num a => a -> a -> a
+Int
1))
                    | Int
j <- [Int
0..Int
nforall a. Num a => a -> a -> a
-Int
1], Int
k <- [Int
0..Int
nforall a. Num a => a -> a -> a
-Int
1], Int
l <- [Int
0..Int
nforall a. Num a => a -> a -> a
-Int
1]]
          tets :: ParamCube -> [Tet]
tets (PC (R, R, R)
c000 (R, R, R)
c001 (R, R, R)
c010 (R, R, R)
c011 (R, R, R)
c100 (R, R, R)
c101 (R, R, R)
c110 (R, R, R)
c111)
              = [Position -> Position -> Position -> Position -> Tet
Tet ((R, R, R) -> Position
g (R, R, R)
c000) ((R, R, R) -> Position
g (R, R, R)
c100) ((R, R, R) -> Position
g (R, R, R)
c010) ((R, R, R) -> Position
g (R, R, R)
c001)
                ,Position -> Position -> Position -> Position -> Tet
Tet ((R, R, R) -> Position
g (R, R, R)
c011) ((R, R, R) -> Position
g (R, R, R)
c111) ((R, R, R) -> Position
g (R, R, R)
c001) ((R, R, R) -> Position
g (R, R, R)
c010)
                ,Position -> Position -> Position -> Position -> Tet
Tet ((R, R, R) -> Position
g (R, R, R)
c110) ((R, R, R) -> Position
g (R, R, R)
c010) ((R, R, R) -> Position
g (R, R, R)
c100) ((R, R, R) -> Position
g (R, R, R)
c111)
                ,Position -> Position -> Position -> Position -> Tet
Tet ((R, R, R) -> Position
g (R, R, R)
c101) ((R, R, R) -> Position
g (R, R, R)
c001) ((R, R, R) -> Position
g (R, R, R)
c111) ((R, R, R) -> Position
g (R, R, R)
c100)
                ,Position -> Position -> Position -> Position -> Tet
Tet ((R, R, R) -> Position
g (R, R, R)
c111) ((R, R, R) -> Position
g (R, R, R)
c100) ((R, R, R) -> Position
g (R, R, R)
c010) ((R, R, R) -> Position
g (R, R, R)
c001)
                ]
      in forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap ParamCube -> [Tet]
tets [ParamCube]
stCubes

type Field a = Position -> a

class AbstractVector a where
    zeroVector :: a
    add   :: a -> a -> a
    scale :: R -> a -> a

sumG :: AbstractVector a => [a] -> a
sumG :: forall a. AbstractVector a => [a] -> a
sumG = forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr forall a. AbstractVector a => a -> a -> a
add forall a. AbstractVector a => a
zeroVector

generalLineIntegral
    :: AbstractVector a => CurveApprox -> Field a -> Curve -> a
generalLineIntegral :: forall a. AbstractVector a => CurveApprox -> Field a -> Curve -> a
generalLineIntegral CurveApprox
approx Field a
f Curve
c
    = forall a. AbstractVector a => [a] -> a
sumG [forall a. AbstractVector a => R -> a -> a
scale (Displacement -> R
magnitude Displacement
dl') (Field a
f Position
r') | (Position
r',Displacement
dl') <- CurveApprox
approx Curve
c]

dottedSurfaceIntegral :: SurfaceApprox -> VectorField -> Surface -> R
dottedSurfaceIntegral :: SurfaceApprox -> VectorField -> Surface -> R
dottedSurfaceIntegral SurfaceApprox
approx VectorField
vF Surface
s
    = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [VectorField
vF Position
r' Displacement -> Displacement -> R
<.> Displacement
da' | (Position
r',Displacement
da') <- SurfaceApprox
approx Surface
s]

electricFluxFromField :: VectorField -> Surface -> R
electricFluxFromField :: VectorField -> Surface -> R
electricFluxFromField = forall a. HasCallStack => a
undefined

electricFluxFromCharge :: ChargeDistribution -> Surface -> R
electricFluxFromCharge :: ChargeDistribution -> Surface -> R
electricFluxFromCharge ChargeDistribution
dist = forall a. HasCallStack => a
undefined ChargeDistribution
dist

eFieldFromSurfaceChargeP :: SurfaceApprox -> ScalarField -> Surface
                         -> VectorField
eFieldFromSurfaceChargeP :: SurfaceApprox -> ScalarField -> Surface -> VectorField
eFieldFromSurfaceChargeP SurfaceApprox
approx ScalarField
sigma Surface
s Position
r
    = [Displacement] -> Displacement
sumV [R -> Position -> VectorField
eFieldFromPointCharge (ScalarField
sigma Position
r' forall a. Num a => a -> a -> a
* Displacement -> R
magnitude Displacement
da') Position
r' Position
r
                | (Position
r',Displacement
da') <- SurfaceApprox
approx Surface
s]

surfaceArea :: Surface -> R
surfaceArea :: Surface -> R
surfaceArea = forall a. HasCallStack => a
undefined

dottedLineIntegral :: CurveApprox -> VectorField -> Curve -> R
dottedLineIntegral :: CurveApprox -> VectorField -> Curve -> R
dottedLineIntegral CurveApprox
approx VectorField
f Curve
c = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [VectorField
f Position
r' Displacement -> Displacement -> R
<.> Displacement
dl' | (Position
r',Displacement
dl') <- CurveApprox
approx Curve
c]

electricPotentialFromField :: VectorField  -- electric field
                           -> ScalarField  -- electric potential
electricPotentialFromField :: VectorField -> ScalarField
electricPotentialFromField VectorField
ef Position
r = forall a. HasCallStack => a
undefined VectorField
ef Position
r