{-# LANGUAGE TypeFamilies #-}
{-# OPTIONS_GHC -Wall #-}
{-# LANGUAGE Safe #-}
module Physics.Learn.Volume
(
Volume(..)
, unitBall
, unitBallCartesian
, centeredBall
, ball
, northernHalfBall
, centeredCylinder
, shiftVolume
, volumeIntegral
)
where
import Data.VectorSpace
( VectorSpace
, InnerSpace
, Scalar
)
import Physics.Learn.CarrotVec
( Vec
, vec
, sumV
, (^+^)
, (^-^)
, (^*)
, (^/)
, (<.>)
, (><)
, magnitude
)
import Physics.Learn.Position
( Position
, Displacement
, Field
, cartesian
, cylindrical
, spherical
, shiftPosition
, displacement
)
data Volume = Volume { Volume -> (Double, Double, Double) -> Position
volumeFunc :: (Double,Double,Double) -> Position
, Volume -> Double
loLimit :: Double
, Volume -> Double
upLimit :: Double
, Volume -> Double -> Double
loCurve :: Double -> Double
, Volume -> Double -> Double
upCurve :: Double -> Double
, Volume -> Double -> Double -> Double
loSurf :: Double -> Double -> Double
, Volume -> Double -> Double -> Double
upSurf :: Double -> Double -> Double
}
unitBall :: Volume
unitBall :: Volume
unitBall = ((Double, Double, Double) -> Position)
-> Double
-> Double
-> (Double -> Double)
-> (Double -> Double)
-> (Double -> Double -> Double)
-> (Double -> Double -> Double)
-> Volume
Volume (Double, Double, Double) -> Position
spherical Double
0 Double
1 (Double -> Double -> Double
forall a b. a -> b -> a
const Double
0) (Double -> Double -> Double
forall a b. a -> b -> a
const Double
forall a. Floating a => a
pi) (\Double
_ Double
_ -> Double
0) (\Double
_ Double
_ -> Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
forall a. Floating a => a
pi)
unitBallCartesian :: Volume
unitBallCartesian :: Volume
unitBallCartesian = ((Double, Double, Double) -> Position)
-> Double
-> Double
-> (Double -> Double)
-> (Double -> Double)
-> (Double -> Double -> Double)
-> (Double -> Double -> Double)
-> Volume
Volume (Double, Double, Double) -> Position
cartesian (-Double
1) Double
1 (\Double
x -> -Double -> Double
sqrtTol (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
x)) (\Double
x -> Double -> Double
sqrtTol (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
x))
(\Double
x Double
y -> -Double -> Double
sqrtTol (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
yDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
y)) (\Double
x Double
y -> Double -> Double
sqrtTol (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
yDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
y))
centeredBall :: Double -> Volume
centeredBall :: Double -> Volume
centeredBall Double
a = ((Double, Double, Double) -> Position)
-> Double
-> Double
-> (Double -> Double)
-> (Double -> Double)
-> (Double -> Double -> Double)
-> (Double -> Double -> Double)
-> Volume
Volume (Double, Double, Double) -> Position
spherical Double
0 Double
a (Double -> Double -> Double
forall a b. a -> b -> a
const Double
0) (Double -> Double -> Double
forall a b. a -> b -> a
const Double
forall a. Floating a => a
pi) (\Double
_ Double
_ -> Double
0) (\Double
_ Double
_ -> Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
forall a. Floating a => a
pi)
ball :: Double
-> Position
-> Volume
ball :: Double -> Position -> Volume
ball Double
a Position
c = ((Double, Double, Double) -> Position)
-> Double
-> Double
-> (Double -> Double)
-> (Double -> Double)
-> (Double -> Double -> Double)
-> (Double -> Double -> Double)
-> Volume
Volume (\(Double
r,Double
th,Double
phi) -> Displacement -> Position -> Position
shiftPosition (Double -> Double -> Double -> Displacement
vec (Double
r Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin Double
th Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos Double
phi) (Double
r Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin Double
th Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin Double
phi) (Double
r Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos Double
th)) Position
c)
Double
0 Double
a (Double -> Double -> Double
forall a b. a -> b -> a
const Double
0) (Double -> Double -> Double
forall a b. a -> b -> a
const Double
forall a. Floating a => a
pi) (\Double
_ Double
_ -> Double
0) (\Double
_ Double
_ -> Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
forall a. Floating a => a
pi)
northernHalfBall :: Volume
northernHalfBall :: Volume
northernHalfBall = ((Double, Double, Double) -> Position)
-> Double
-> Double
-> (Double -> Double)
-> (Double -> Double)
-> (Double -> Double -> Double)
-> (Double -> Double -> Double)
-> Volume
Volume (Double, Double, Double) -> Position
spherical Double
0 Double
1 (Double -> Double -> Double
forall a b. a -> b -> a
const Double
0) (Double -> Double -> Double
forall a b. a -> b -> a
const (Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2)) (\Double
_ Double
_ -> Double
0) (\Double
_ Double
_ -> Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
forall a. Floating a => a
pi)
centeredCylinder :: Double
-> Double
-> Volume
centeredCylinder :: Double -> Double -> Volume
centeredCylinder Double
r Double
h = ((Double, Double, Double) -> Position)
-> Double
-> Double
-> (Double -> Double)
-> (Double -> Double)
-> (Double -> Double -> Double)
-> (Double -> Double -> Double)
-> Volume
Volume (Double, Double, Double) -> Position
cylindrical Double
0 Double
r (Double -> Double -> Double
forall a b. a -> b -> a
const Double
0) (Double -> Double -> Double
forall a b. a -> b -> a
const (Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
forall a. Floating a => a
pi)) (\Double
_ Double
_ -> Double
0) (\Double
_ Double
_ -> Double
h)
volumeIntegral :: (VectorSpace v, Scalar v ~ Double) =>
Int
-> Int
-> Int
-> Field v
-> Volume
-> v
volumeIntegral :: forall v.
(VectorSpace v, Scalar v ~ Double) =>
Int -> Int -> Int -> Field v -> Volume -> v
volumeIntegral Int
n1 Int
n2 Int
n3 Field v
field (Volume (Double, Double, Double) -> Position
f Double
s_l Double
s_u Double -> Double
t_l Double -> Double
t_u Double -> Double -> Double
u_l Double -> Double -> Double
u_u)
= [v] -> v
forall (f :: * -> *) v. (Foldable f, AdditiveGroup v) => f v -> v
sumV ([v] -> v) -> [v] -> v
forall a b. (a -> b) -> a -> b
$ ([v] -> v) -> [[v]] -> [v]
forall a b. (a -> b) -> [a] -> [b]
map [v] -> v
forall (f :: * -> *) v. (Foldable f, AdditiveGroup v) => f v -> v
sumV ([[v]] -> [v]) -> [[v]] -> [v]
forall a b. (a -> b) -> a -> b
$ ([[v]] -> [v]) -> [[[v]]] -> [[v]]
forall a b. (a -> b) -> [a] -> [b]
map (([v] -> v) -> [[v]] -> [v]
forall a b. (a -> b) -> [a] -> [b]
map [v] -> v
forall (f :: * -> *) v. (Foldable f, AdditiveGroup v) => f v -> v
sumV) ((v -> Double -> v) -> [[[v]]] -> [[[Double]]] -> [[[v]]]
forall a b c. (a -> b -> c) -> [[[a]]] -> [[[b]]] -> [[[c]]]
zipCubeWith v -> Double -> v
forall v s. (VectorSpace v, s ~ Scalar v) => v -> s -> v
(^*) [[[v]]]
aveVals [[[Double]]]
volumes)
where
pts :: [[[Position]]]
pts = [[[(Double, Double, Double) -> Position
f (Double
s,Double
t,Double
u) | Double
u <- Int -> Double -> Double -> [Double]
linSpaced Int
n3 (Double -> Double -> Double
u_l Double
s Double
t) (Double -> Double -> Double
u_u Double
s Double
t) ] | Double
t <- Int -> Double -> Double -> [Double]
linSpaced Int
n2 (Double -> Double
t_l Double
s) (Double -> Double
t_u Double
s)] | Double
s <- Int -> Double -> Double -> [Double]
linSpaced Int
n1 Double
s_l Double
s_u]
volumes :: [[[Double]]]
volumes = ([[Displacement]]
-> [[Displacement]] -> [[Displacement]] -> [[Double]])
-> [[[Displacement]]]
-> [[[Displacement]]]
-> [[[Displacement]]]
-> [[[Double]]]
forall a b c d. (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
zipWith3 (([Displacement] -> [Displacement] -> [Displacement] -> [Double])
-> [[Displacement]]
-> [[Displacement]]
-> [[Displacement]]
-> [[Double]]
forall a b c d. (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
zipWith3 ((Displacement -> Displacement -> Displacement -> Double)
-> [Displacement] -> [Displacement] -> [Displacement] -> [Double]
forall a b c d. (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
zipWith3 (\Displacement
du Displacement
dv Displacement
dw -> Displacement
du Displacement -> Displacement -> Scalar Displacement
forall v. InnerSpace v => v -> v -> Scalar v
<.> (Displacement
dv Displacement -> Displacement -> Displacement
>< Displacement
dw)))) [[[Displacement]]]
dus [[[Displacement]]]
dvs [[[Displacement]]]
dws
dus :: [[[Displacement]]]
dus = ([[[Position]]] -> [[[Position]]] -> [[[Displacement]]])
-> ([[[Position]]], [[[Position]]]) -> [[[Displacement]]]
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry [[[Position]]] -> [[[Position]]] -> [[[Displacement]]]
zipSub3 ([[[Position]]] -> ([[[Position]]], [[[Position]]])
forall a. [a] -> ([a], [a])
shift1 [[[Position]]]
pts)
dvs :: [[[Displacement]]]
dvs = ([[[Position]]] -> [[[Position]]] -> [[[Displacement]]])
-> ([[[Position]]], [[[Position]]]) -> [[[Displacement]]]
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry [[[Position]]] -> [[[Position]]] -> [[[Displacement]]]
zipSub3 ([[[Position]]] -> ([[[Position]]], [[[Position]]])
forall a. [[a]] -> ([[a]], [[a]])
shift2 [[[Position]]]
pts)
dws :: [[[Displacement]]]
dws = ([[[Position]]] -> [[[Position]]] -> [[[Displacement]]])
-> ([[[Position]]], [[[Position]]]) -> [[[Displacement]]]
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry [[[Position]]] -> [[[Position]]] -> [[[Displacement]]]
zipSub3 ([[[Position]]] -> ([[[Position]]], [[[Position]]])
forall a. [[[a]]] -> ([[[a]]], [[[a]]])
shift3 [[[Position]]]
pts)
vals :: [[[v]]]
vals = ([[Position]] -> [[v]]) -> [[[Position]]] -> [[[v]]]
forall a b. (a -> b) -> [a] -> [b]
map (([Position] -> [v]) -> [[Position]] -> [[v]]
forall a b. (a -> b) -> [a] -> [b]
map (Field v -> [Position] -> [v]
forall a b. (a -> b) -> [a] -> [b]
map Field v
field)) [[[Position]]]
pts
aveVals :: [[[v]]]
aveVals = ((([[[v]]] -> [[[v]]] -> [[[v]]]) -> ([[[v]]], [[[v]]]) -> [[[v]]]
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry [[[v]]] -> [[[v]]] -> [[[v]]]
forall v.
(VectorSpace v, Scalar v ~ Double) =>
[[[v]]] -> [[[v]]] -> [[[v]]]
zipAve3 (([[[v]]], [[[v]]]) -> [[[v]]])
-> ([[[v]]] -> ([[[v]]], [[[v]]])) -> [[[v]]] -> [[[v]]]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [[[v]]] -> ([[[v]]], [[[v]]])
forall a. [a] -> ([a], [a])
shift1) ([[[v]]] -> [[[v]]]) -> ([[[v]]] -> [[[v]]]) -> [[[v]]] -> [[[v]]]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (([[[v]]] -> [[[v]]] -> [[[v]]]) -> ([[[v]]], [[[v]]]) -> [[[v]]]
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry [[[v]]] -> [[[v]]] -> [[[v]]]
forall v.
(VectorSpace v, Scalar v ~ Double) =>
[[[v]]] -> [[[v]]] -> [[[v]]]
zipAve3 (([[[v]]], [[[v]]]) -> [[[v]]])
-> ([[[v]]] -> ([[[v]]], [[[v]]])) -> [[[v]]] -> [[[v]]]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [[[v]]] -> ([[[v]]], [[[v]]])
forall a. [[a]] -> ([[a]], [[a]])
shift2) ([[[v]]] -> [[[v]]]) -> ([[[v]]] -> [[[v]]]) -> [[[v]]] -> [[[v]]]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (([[[v]]] -> [[[v]]] -> [[[v]]]) -> ([[[v]]], [[[v]]]) -> [[[v]]]
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry [[[v]]] -> [[[v]]] -> [[[v]]]
forall v.
(VectorSpace v, Scalar v ~ Double) =>
[[[v]]] -> [[[v]]] -> [[[v]]]
zipAve3 (([[[v]]], [[[v]]]) -> [[[v]]])
-> ([[[v]]] -> ([[[v]]], [[[v]]])) -> [[[v]]] -> [[[v]]]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [[[v]]] -> ([[[v]]], [[[v]]])
forall a. [[[a]]] -> ([[[a]]], [[[a]]])
shift3)) [[[v]]]
vals
zipCubeWith :: (a -> b -> c) -> [[[a]]] -> [[[b]]] -> [[[c]]]
zipCubeWith :: forall a b c. (a -> b -> c) -> [[[a]]] -> [[[b]]] -> [[[c]]]
zipCubeWith = ([[a]] -> [[b]] -> [[c]]) -> [[[a]]] -> [[[b]]] -> [[[c]]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (([[a]] -> [[b]] -> [[c]]) -> [[[a]]] -> [[[b]]] -> [[[c]]])
-> ((a -> b -> c) -> [[a]] -> [[b]] -> [[c]])
-> (a -> b -> c)
-> [[[a]]]
-> [[[b]]]
-> [[[c]]]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ([a] -> [b] -> [c]) -> [[a]] -> [[b]] -> [[c]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (([a] -> [b] -> [c]) -> [[a]] -> [[b]] -> [[c]])
-> ((a -> b -> c) -> [a] -> [b] -> [c])
-> (a -> b -> c)
-> [[a]]
-> [[b]]
-> [[c]]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a -> b -> c) -> [a] -> [b] -> [c]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith
zipSub3 :: [[[Position]]] -> [[[Position]]] -> [[[Vec]]]
zipSub3 :: [[[Position]]] -> [[[Position]]] -> [[[Displacement]]]
zipSub3 = (Position -> Position -> Displacement)
-> [[[Position]]] -> [[[Position]]] -> [[[Displacement]]]
forall a b c. (a -> b -> c) -> [[[a]]] -> [[[b]]] -> [[[c]]]
zipCubeWith Position -> Position -> Displacement
displacement
zipAve3 :: (VectorSpace v, Scalar v ~ Double) => [[[v]]] -> [[[v]]] -> [[[v]]]
zipAve3 :: forall v.
(VectorSpace v, Scalar v ~ Double) =>
[[[v]]] -> [[[v]]] -> [[[v]]]
zipAve3 = (v -> v -> v) -> [[[v]]] -> [[[v]]] -> [[[v]]]
forall a b c. (a -> b -> c) -> [[[a]]] -> [[[b]]] -> [[[c]]]
zipCubeWith v -> v -> v
forall v. (VectorSpace v, Scalar v ~ Double) => v -> v -> v
ave
shift1 :: [a] -> ([a],[a])
shift1 :: forall a. [a] -> ([a], [a])
shift1 [a]
pts = ([a]
pts, [a] -> [a]
forall a. HasCallStack => [a] -> [a]
tail [a]
pts)
shift2 :: [[a]] -> ([[a]],[[a]])
shift2 :: forall a. [[a]] -> ([[a]], [[a]])
shift2 [[a]]
pts2d = ([[a]]
pts2d, ([a] -> [a]) -> [[a]] -> [[a]]
forall a b. (a -> b) -> [a] -> [b]
map [a] -> [a]
forall a. HasCallStack => [a] -> [a]
tail [[a]]
pts2d)
shift3 :: [[[a]]] -> ([[[a]]],[[[a]]])
shift3 :: forall a. [[[a]]] -> ([[[a]]], [[[a]]])
shift3 [[[a]]]
pts3d = ([[[a]]]
pts3d, ([[a]] -> [[a]]) -> [[[a]]] -> [[[a]]]
forall a b. (a -> b) -> [a] -> [b]
map (([a] -> [a]) -> [[a]] -> [[a]]
forall a b. (a -> b) -> [a] -> [b]
map [a] -> [a]
forall a. HasCallStack => [a] -> [a]
tail) [[[a]]]
pts3d)
linSpaced :: Int -> Double -> Double -> [Double]
linSpaced :: Int -> Double -> Double -> [Double]
linSpaced Int
n Double
a Double
b
| Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
b = let dx :: Double
dx = (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
in [Double
a,Double
aDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
dx..Double
b]
| Double
a Double -> Double -> Bool
forall v. (InnerSpace v, Scalar v ~ Double) => v -> v -> Bool
~~ Double
b = [Double -> Double -> Double
forall v. (VectorSpace v, Scalar v ~ Double) => v -> v -> v
ave Double
a Double
b]
| Bool
otherwise = [Char] -> [Double]
forall a. HasCallStack => [Char] -> a
error ([Char] -> [Double]) -> [Char] -> [Double]
forall a b. (a -> b) -> a -> b
$ [Char]
"linSpaced: lower limit greater than upper limit: (n,a,b) = " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ (Int, Double, Double) -> [Char]
forall a. Show a => a -> [Char]
show (Int
n,Double
a,Double
b)
(~~) :: (InnerSpace v, Scalar v ~ Double) => v -> v -> Bool
v
a ~~ :: forall v. (InnerSpace v, Scalar v ~ Double) => v -> v -> Bool
~~ v
b = v -> Double
forall v s. (InnerSpace v, s ~ Scalar v, Floating s) => v -> s
magnitude (v
a v -> v -> v
forall v. AdditiveGroup v => v -> v -> v
^-^ v
b) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
tolerance
tolerance :: Double
tolerance :: Double
tolerance = Double
1e-10
ave :: (VectorSpace v, Scalar v ~ Double) => v -> v -> v
ave :: forall v. (VectorSpace v, Scalar v ~ Double) => v -> v -> v
ave v
v1 v
v2 = (v
v1 v -> v -> v
forall v. AdditiveGroup v => v -> v -> v
^+^ v
v2) v -> Double -> v
forall v s.
(VectorSpace v, s ~ Scalar v, Fractional s) =>
v -> s -> v
^/ Double
2
sqrtTol :: Double -> Double
sqrtTol :: Double -> Double
sqrtTol Double
x
| Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0 = Double -> Double
forall a. Floating a => a -> a
sqrt Double
x
| Double -> Double
forall a. Num a => a -> a
abs Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
tolerance = Double
0
| Bool
otherwise = [Char] -> Double
forall a. HasCallStack => [Char] -> a
error ([Char]
"sqrtTol: I can't take the sqrt of " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ Double -> [Char]
forall a. Show a => a -> [Char]
show Double
x)
shiftVolume :: Displacement -> Volume -> Volume
shiftVolume :: Displacement -> Volume -> Volume
shiftVolume Displacement
d (Volume (Double, Double, Double) -> Position
f Double
sl Double
su Double -> Double
tl Double -> Double
tu Double -> Double -> Double
ul Double -> Double -> Double
uu)
= ((Double, Double, Double) -> Position)
-> Double
-> Double
-> (Double -> Double)
-> (Double -> Double)
-> (Double -> Double -> Double)
-> (Double -> Double -> Double)
-> Volume
Volume (Displacement -> Position -> Position
shiftPosition Displacement
d (Position -> Position)
-> ((Double, Double, Double) -> Position)
-> (Double, Double, Double)
-> Position
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Double, Double, Double) -> Position
f) Double
sl Double
su Double -> Double
tl Double -> Double
tu Double -> Double -> Double
ul Double -> Double -> Double
uu