{-# LANGUAGE TypeFamilies #-}
{-# OPTIONS_GHC -Wall #-}
{-# LANGUAGE Safe #-}

{- | 
Module      :  Physics.Learn.Volume
Copyright   :  (c) Scott N. Walck 2012-2019
License     :  BSD3 (see LICENSE)
Maintainer  :  Scott N. Walck <walck@lvc.edu>
Stability   :  experimental

This module contains functions for working with 'Volume's
and volume integrals over 'Volume's.
-}

module Physics.Learn.Volume
    (
    -- * Volumes
      Volume(..)
    , unitBall
    , unitBallCartesian
    , centeredBall
    , ball
    , northernHalfBall
    , centeredCylinder
    , shiftVolume
    -- * Volume Integral
    , 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
    )

-- | Volume is a parametrized function from three parameters to space,
--   lower and upper limits on the first parameter,
--   lower and upper limits for the second parameter
--   (expressed as functions of the first parameter),
--   and lower and upper limits for the third parameter
--   (expressed as functions of the first and second parameters).
data Volume = Volume { Volume -> (Double, Double, Double) -> Position
volumeFunc :: (Double,Double,Double) -> Position  -- ^ function from 3 parameters to space
                     , Volume -> Double
loLimit    :: Double                      -- ^ s_a
                     , Volume -> Double
upLimit    :: Double                      -- ^ s_b
                     , Volume -> Double -> Double
loCurve    :: Double -> Double            -- ^ t_a(s)
                     , Volume -> Double -> Double
upCurve    :: Double -> Double            -- ^ t_b(s)
                     , Volume -> Double -> Double -> Double
loSurf     :: Double -> Double -> Double  -- ^ u_a(s,t)
                     , Volume -> Double -> Double -> Double
upSurf     :: Double -> Double -> Double  -- ^ u_b(s,t)
                     }

-- | A unit ball, centered at the origin.
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 (forall a b. a -> b -> a
const Double
0) (forall a b. a -> b -> a
const forall a. Floating a => a
pi) (\Double
_ Double
_ -> Double
0) (\Double
_ Double
_ -> Double
2forall a. Num a => a -> a -> a
*forall a. Floating a => a
pi)

-- | A unit ball, centered at the origin.  Specified in Cartesian coordinates.
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 forall a. Num a => a -> a -> a
- Double
xforall a. Num a => a -> a -> a
*Double
x)) (\Double
x -> Double -> Double
sqrtTol (Double
1 forall a. Num a => a -> a -> a
- Double
xforall a. Num a => a -> a -> a
*Double
x))
                    (\Double
x Double
y -> -Double -> Double
sqrtTol (Double
1 forall a. Num a => a -> a -> a
- Double
xforall a. Num a => a -> a -> a
*Double
x forall a. Num a => a -> a -> a
- Double
yforall a. Num a => a -> a -> a
*Double
y)) (\Double
x Double
y -> Double -> Double
sqrtTol (Double
1 forall a. Num a => a -> a -> a
- Double
xforall a. Num a => a -> a -> a
*Double
x forall a. Num a => a -> a -> a
- Double
yforall a. Num a => a -> a -> a
*Double
y))

-- | A ball with given radius, centered at the origin.
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 (forall a b. a -> b -> a
const Double
0) (forall a b. a -> b -> a
const forall a. Floating a => a
pi) (\Double
_ Double
_ -> Double
0) (\Double
_ Double
_ -> Double
2forall a. Num a => a -> a -> a
*forall a. Floating a => a
pi)

-- | Ball with given radius and center.
ball :: Double    -- ^ radius
     -> Position  -- ^ center
     -> Volume    -- ^ ball with given radius and center
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 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sin Double
th forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
cos Double
phi) (Double
r forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sin Double
th forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sin Double
phi) (Double
r forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
cos Double
th)) Position
c)
           Double
0 Double
a (forall a b. a -> b -> a
const Double
0) (forall a b. a -> b -> a
const forall a. Floating a => a
pi) (\Double
_ Double
_ -> Double
0) (\Double
_ Double
_ -> Double
2forall a. Num a => a -> a -> a
*forall a. Floating a => a
pi)

-- | Upper half ball, unit radius, centered at origin.
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 (forall a b. a -> b -> a
const Double
0) (forall a b. a -> b -> a
const (forall a. Floating a => a
piforall a. Fractional a => a -> a -> a
/Double
2)) (\Double
_ Double
_ -> Double
0) (\Double
_ Double
_ -> Double
2forall a. Num a => a -> a -> a
*forall a. Floating a => a
pi)

-- | Cylinder with given radius and height.  Circular base of the cylinder
--   is centered at the origin.  Circular top of the cylinder lies in plane z = h.
centeredCylinder :: Double  -- radius
                 -> Double  -- height
                 -> Volume  -- cylinder
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 (forall a b. a -> b -> a
const Double
0) (forall a b. a -> b -> a
const (Double
2forall a. Num a => a -> a -> a
*forall a. Floating a => a
pi)) (\Double
_ Double
_ -> Double
0) (\Double
_ Double
_ -> Double
h)

-- | A volume integral
volumeIntegral :: (VectorSpace v, Scalar v ~ Double) =>
                  Int          -- ^ number of intervals for first parameter   (s)
               -> Int          -- ^ number of intervals for second parameter  (t)
               -> Int          -- ^ number of intervals for third parameter   (u)
               -> Field v      -- ^ scalar or vector field
               -> Volume       -- ^ the volume
               -> v            -- ^ scalar or vector result
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)
    = forall (f :: * -> *) v. (Foldable f, AdditiveGroup v) => f v -> v
sumV forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map forall (f :: * -> *) v. (Foldable f, AdditiveGroup v) => f v -> v
sumV forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map (forall a b. (a -> b) -> [a] -> [b]
map forall (f :: * -> *) v. (Foldable f, AdditiveGroup v) => f v -> v
sumV) (forall a b c. (a -> b -> c) -> [[[a]]] -> [[[b]]] -> [[[c]]]
zipCubeWith 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 = forall a b c d. (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
zipWith3 (forall a b c d. (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
zipWith3 (forall a b c d. (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
zipWith3 (\Displacement
du Displacement
dv Displacement
dw -> Displacement
du forall v. InnerSpace v => v -> v -> Scalar v
<.> (Displacement
dv Displacement -> Displacement -> Displacement
>< Displacement
dw)))) [[[Displacement]]]
dus [[[Displacement]]]
dvs [[[Displacement]]]
dws
        dus :: [[[Displacement]]]
dus = forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry [[[Position]]] -> [[[Position]]] -> [[[Displacement]]]
zipSub3 (forall a. [a] -> ([a], [a])
shift1 [[[Position]]]
pts)
        dvs :: [[[Displacement]]]
dvs = forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry [[[Position]]] -> [[[Position]]] -> [[[Displacement]]]
zipSub3 (forall a. [[a]] -> ([[a]], [[a]])
shift2 [[[Position]]]
pts)
        dws :: [[[Displacement]]]
dws = forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry [[[Position]]] -> [[[Position]]] -> [[[Displacement]]]
zipSub3 (forall a. [[[a]]] -> ([[[a]]], [[[a]]])
shift3 [[[Position]]]
pts)
        vals :: [[[v]]]
vals = forall a b. (a -> b) -> [a] -> [b]
map (forall a b. (a -> b) -> [a] -> [b]
map (forall a b. (a -> b) -> [a] -> [b]
map Field v
field)) [[[Position]]]
pts
        aveVals :: [[[v]]]
aveVals = ((forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry forall v.
(VectorSpace v, Scalar v ~ Double) =>
[[[v]]] -> [[[v]]] -> [[[v]]]
zipAve3 forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. [a] -> ([a], [a])
shift1) forall b c a. (b -> c) -> (a -> b) -> a -> c
. (forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry forall v.
(VectorSpace v, Scalar v ~ Double) =>
[[[v]]] -> [[[v]]] -> [[[v]]]
zipAve3 forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. [[a]] -> ([[a]], [[a]])
shift2) forall b c a. (b -> c) -> (a -> b) -> a -> c
. (forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry forall v.
(VectorSpace v, Scalar v ~ Double) =>
[[[v]]] -> [[[v]]] -> [[[v]]]
zipAve3 forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. [[[a]]] -> ([[[a]]], [[[a]]])
shift3)) [[[v]]]
vals

-- zipSquareWith :: (a -> b -> c) -> [[a]] -> [[b]] -> [[c]]
-- zipSquareWith = zipWith . zipWith

zipCubeWith :: (a -> b -> c) -> [[[a]]] -> [[[b]]] -> [[[c]]]
zipCubeWith :: forall a b c. (a -> b -> c) -> [[[a]]] -> [[[b]]] -> [[[c]]]
zipCubeWith = forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith

-- zipSub :: [Vec] -> [Vec] -> [Vec]
-- zipSub = zipWith (^-^)

-- zipSub2 :: [[Vec]] -> [[Vec]] -> [[Vec]]
-- zipSub2 = zipWith $ zipWith (^-^)

zipSub3 :: [[[Position]]] -> [[[Position]]] -> [[[Vec]]]
zipSub3 :: [[[Position]]] -> [[[Position]]] -> [[[Displacement]]]
zipSub3 = 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 = forall a b c. (a -> b -> c) -> [[[a]]] -> [[[b]]] -> [[[c]]]
zipCubeWith 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, forall a. [a] -> [a]
tail [a]
pts)

shift2 :: [[a]] -> ([[a]],[[a]])
shift2 :: forall a. [[a]] -> ([[a]], [[a]])
shift2 [[a]]
pts2d = ([[a]]
pts2d, forall a b. (a -> b) -> [a] -> [b]
map forall a. [a] -> [a]
tail [[a]]
pts2d)

shift3 :: [[[a]]] -> ([[[a]]],[[[a]]])
shift3 :: forall a. [[[a]]] -> ([[[a]]], [[[a]]])
shift3 [[[a]]]
pts3d = ([[[a]]]
pts3d, forall a b. (a -> b) -> [a] -> [b]
map (forall a b. (a -> b) -> [a] -> [b]
map forall a. [a] -> [a]
tail) [[[a]]]
pts3d)

-- | n+1 points
linSpaced :: Int -> Double -> Double -> [Double]
linSpaced :: Int -> Double -> Double -> [Double]
linSpaced Int
n Double
a Double
b
    | Double
a forall a. Ord a => a -> a -> Bool
< Double
b      = let dx :: Double
dx = (Double
b forall a. Num a => a -> a -> a
- Double
a) forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
                   in [Double
a,Double
aforall a. Num a => a -> a -> a
+Double
dx..Double
b]
    | Double
a forall v. (InnerSpace v, Scalar v ~ Double) => v -> v -> Bool
~~ Double
b     = [forall v. (VectorSpace v, Scalar v ~ Double) => v -> v -> v
ave Double
a Double
b]
    | Bool
otherwise  = forall a. HasCallStack => [Char] -> a
error forall a b. (a -> b) -> a -> b
$ [Char]
"linSpaced:  lower limit greater than upper limit:  (n,a,b) = " forall a. [a] -> [a] -> [a]
++ 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 = forall v s. (InnerSpace v, s ~ Scalar v, Floating s) => v -> s
magnitude (v
a forall v. AdditiveGroup v => v -> v -> v
^-^ v
b) 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 forall v. AdditiveGroup v => v -> v -> v
^+^ v
v2) 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 forall a. Ord a => a -> a -> Bool
>= Double
0              = forall a. Floating a => a -> a
sqrt Double
x
    | forall a. Num a => a -> a
abs Double
x forall a. Ord a => a -> a -> Bool
<= Double
tolerance  = Double
0
    | Bool
otherwise           = forall a. HasCallStack => [Char] -> a
error ([Char]
"sqrtTol:  I can't take the sqrt of " forall a. [a] -> [a] -> [a]
++ forall a. Show a => a -> [Char]
show Double
x)

-- | Shift a volume by a displacement.
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 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