module Numerics.VectorAnalysis (

    -- * Operators with default settings
      grad
    , divergence
    , curl
    , curlZ
    , laplace

    -- * Operators with configurable step width
    , gradH
    , divH
    , curlH
    , curlZH
    , laplaceH
) where

import Geometry.Core

-- | A good standard value to use as a step size for taking derivatives.
standardH :: Double
standardH :: Double
standardH = Double
1e-3

-- | Gradient with predefined sampling distance.
--
-- \[
-- \text{grad}(f) = \begin{pmatrix}
--         \partial_x f \\
--         \partial_y f
--     \end{pmatrix}
-- \]
--
-- The fire-and-forget version of 'gradH'.
grad :: (Vec2 -> Double) -> Vec2 -> Vec2
grad :: (Vec2 -> Double) -> Vec2 -> Vec2
grad = Double -> (Vec2 -> Double) -> Vec2 -> Vec2
gradH Double
standardH

-- | Divergence with predefined sampling distance. Named to avoid clashing with
-- integer 'div'ision.
--
-- \[
-- \text{div}(f) = \partial_xf_x + \partial_y f_y
-- \]
--
-- The fire-and-forget version of 'divH'.
divergence :: (Vec2 -> Vec2) -> Vec2 -> Double
divergence :: (Vec2 -> Vec2) -> Vec2 -> Double
divergence = Double -> (Vec2 -> Vec2) -> Vec2 -> Double
divH Double
standardH

-- | Curl with predefined sampling distance. Note that in two dimensions, the curl
-- is simply a number. For a different 2D adaptation of the \(\text{curl}\)
-- operator, see 'curlZ'.
--
-- \[
-- \text{curl}(f) = \partial_x f_y - \partial_y f_x
-- \]
--
-- The fire-and-forget version of 'curlH'.
curl :: (Vec2 -> Vec2) -> Vec2 -> Double
curl :: (Vec2 -> Vec2) -> Vec2 -> Double
curl = Double -> (Vec2 -> Vec2) -> Vec2 -> Double
curlH Double
standardH

-- | Laplacian with predefined sampling distance.
--
-- \[
-- \nabla^2f = \partial^2_x f + \partial^2_y f
-- \]
--
-- The fire-and-forget version of 'laplaceH'.
laplace :: (Vec2 -> Double) -> Vec2 -> Double
laplace :: (Vec2 -> Double) -> Vec2 -> Double
laplace = Double -> (Vec2 -> Double) -> Vec2 -> Double
laplaceH Double
standardH

-- | Gradient with customizable sampling distance. The configurable version of 'grad'.
gradH
    :: Double -- ^ \(h\) as in \(\frac{f(x+h)-f(x)}h\)
    -> (Vec2 -> Double)
    -> (Vec2 -> Vec2)
gradH :: Double -> (Vec2 -> Double) -> Vec2 -> Vec2
gradH Double
h Vec2 -> Double
f = \Vec2
x ->
    let f_x :: Double
f_x = Vec2 -> Double
f Vec2
x
    in Double -> Double -> Vec2
Vec2 (Vec2 -> Double
f (Vec2
x Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
+. Double -> Double -> Vec2
Vec2 Double
h Double
0) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
f_x) (Vec2 -> Double
f (Vec2
x Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
+. Double -> Double -> Vec2
Vec2 Double
0 Double
h) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
f_x) Vec2 -> Double -> Vec2
forall v. VectorSpace v => v -> Double -> v
/. Double
h

-- | Divergence with customizable sampling distance. The configurable version of 'divergence'.
divH
    :: Double -- ^ \(h\) as in \(\frac{f(x+h)-f(x)}h\)
    -> (Vec2 -> Vec2)
    -> (Vec2 -> Double)
divH :: Double -> (Vec2 -> Vec2) -> Vec2 -> Double
divH Double
h Vec2 -> Vec2
f = \v :: Vec2
v@(Vec2 Double
x Double
y) ->
    let f_v :: Vec2
f_v = Vec2 -> Vec2
f Vec2
v
        Vec2 Double
dx Double
_ = Vec2 -> Vec2
f (Double -> Double -> Vec2
Vec2 (Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
h) Double
y) Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
-. Vec2
f_v
        Vec2 Double
_ Double
dy = Vec2 -> Vec2
f (Double -> Double -> Vec2
Vec2 Double
x (Double
yDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
h)) Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
-. Vec2
f_v
    in (Double
dxDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
dy)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
h

-- | Curl with customizable sampling distance. The configurable version of 'curl'.
curlH
    :: Double -- ^ \(h\) as in \(\frac{f(x+h)-f(x)}h\)
    -> (Vec2 -> Vec2)
    -> (Vec2 -> Double)
curlH :: Double -> (Vec2 -> Vec2) -> Vec2 -> Double
curlH Double
h Vec2 -> Vec2
f = \v :: Vec2
v@(Vec2 Double
x Double
y) ->
    let f_v :: Vec2
f_v = Vec2 -> Vec2
f Vec2
v
        Vec2 Double
dy_fx Double
_ = Vec2 -> Vec2
f (Double -> Double -> Vec2
Vec2 Double
x (Double
yDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
h)) Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
-. Vec2
f_v
        Vec2 Double
_ Double
dx_fy = Vec2 -> Vec2
f (Double -> Double -> Vec2
Vec2 (Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
h) Double
y) Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
-. Vec2
f_v
    in (Double
dx_fyDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
dy_fx) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
h

-- | Laplacian with customizable sampling distance. The configurable version of 'laplace'.
laplaceH
    :: Double -- ^ \(h\) as in \(\frac{f(x+h)-f(x)}h\)
    -> (Vec2 -> Double)
    -> (Vec2 -> Double)
laplaceH :: Double -> (Vec2 -> Double) -> Vec2 -> Double
laplaceH Double
h = Double -> (Vec2 -> Vec2) -> Vec2 -> Double
divH Double
h ((Vec2 -> Vec2) -> Vec2 -> Double)
-> ((Vec2 -> Double) -> Vec2 -> Vec2)
-> (Vec2 -> Double)
-> Vec2
-> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> (Vec2 -> Double) -> Vec2 -> Vec2
gradH Double
h

-- | Curl of a purely-z-component 3D vector field, which is another common way
-- (than 'curl') to implement a two-dimensional version of \(\text{curl}\):
--
-- \[
-- \text{curl}_z(\psi)
--     = \left[\text{curl}_{3D} \begin{pmatrix}0\\0\\\psi\end{pmatrix}\right]_{\text{2D part}}
--     = \begin{pmatrix}\partial_y\psi\\-\partial_x\psi\end{pmatrix}
-- \]
--
-- Because of this, the result is always divergence-free, because \(\forall f. \text{div}(\text{curl}(f))=0\).
--
-- Using 'curlZ' to create 'divergence'-free flow fields is similar to using 'grad'
-- to create 'curl'-free force fields. The argument for 'curlZ' is known as the
-- vector or stream potential in fluid dynamics, with the resulting vector field
-- being known as a
-- [stream function](https://en.wikipedia.org/wiki/Stream_function).
--
-- The fire-and-forget version of 'curlZH'.
curlZ :: (Vec2 -> Double) -> Vec2 -> Vec2
curlZ :: (Vec2 -> Double) -> Vec2 -> Vec2
curlZ Vec2 -> Double
f = Double -> (Vec2 -> Double) -> Vec2 -> Vec2
curlZH Double
standardH Vec2 -> Double
f

-- | Curl of the z component of a 3D vector field, with customizable sampling
-- distance. The configurable version of 'curlZ'.
curlZH
    :: Double -- ^ \(h\) as in \(\frac{f(x+h)-f(x)}h\)
    -> (Vec2 -> Double)
    -> Vec2
    -> Vec2
curlZH :: Double -> (Vec2 -> Double) -> Vec2 -> Vec2
curlZH Double
h Vec2 -> Double
f = \Vec2
x ->
    let f_x :: Double
f_x = Vec2 -> Double
f Vec2
x
    in Double -> Double -> Vec2
Vec2 (Vec2 -> Double
f (Vec2
x Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
+. Double -> Double -> Vec2
Vec2 Double
0 Double
h) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
f_x) (- Vec2 -> Double
f (Vec2
x Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
+. Double -> Double -> Vec2
Vec2 Double
h Double
0) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
f_x) Vec2 -> Double -> Vec2
forall v. VectorSpace v => v -> Double -> v
/. Double
h
        -- NB: this is just 90° rotated grad