module Geometry.Processes.Geodesics (geodesicEquation) where



import Geometry.Core



-- | A geodesic is a path that takes no local detours: each move to a new point is
-- done so that no other move would be faster. The shortest path between two points
-- is always a geodesic.
--
-- This function allows creating the geodesic differential equation, suitable for
-- using in ODE solvers such as
-- 'Numerics.DifferentialEquation.rungeKuttaAdaptiveStep'.
--
-- Using this is very simple, the equation behind it is simple depending on your
-- knowledge of differential geometry, and implementing it is a quite interesting
-- exercise in algorithmic optimization.
--
-- \[
-- \begin{align}
--     \ddot v^i &= \Gamma^i_{kl}\dot v^k\dot v^l \\
--     \Gamma^i_{kl} &= \frac12 g^{im} (g_{mk,l}+g_{ml,k}-g_{kl,m}) \\
--     g(f) &= \begin{pmatrix}1+\partial_{xx}f & \partial_xf\,\partial_yf \\ \partial_xf\,\partial_yf & 1+\partial_{yy}f\end{pmatrix}
-- \end{align}
-- \]
--
-- == Example: hill and valley
--
-- We can use the following variation of the
-- [Cauchy distribution](https://en.wikipedia.org/wiki/Cauchy_distribution) to
-- build ourselves a hill with width parameter \(\gamma\),
--
-- \[
-- \text{Hill}_{\gamma}(\mathbf v_{\text{center}}, \mathbf v)
--     = \left(1+\left\| \frac{\mathbf v-\mathbf v_{\text{center}}}{\gamma}\right\|^2\right)^{-1}
-- \]
--
-- We can then pleace two of those, one at the top right with height 100 and one at
-- the bottom left with height -33.
--
-- This is what a family of trajectories passing through our terrain looks like:
--
-- @
-- (w,h) = … -- Canvas size
-- hill g v0 v = 1 \/ (1 + ('normSquare' (v-.v0) \/ g^2))
-- geometry v = 100 * hill 55 ('Vec2' (2\/3*w) (1\/3*h)) v - 33 * hill 55 ('Vec2' (1\/3*w) (2\/3*h)) v
--
-- startAngle = 'deg' 45
-- solution = 'Numerics.DifferentialEquation.rungeKuttaAdaptiveStep' ('geodesicEquation' geometry) ('zero', 'polar' startAngle 1) t0 dt0 tolNorm tol)
--   where
--     t0 = 0
--     dt0 = 1
--     tolNorm (x,v) = 'max' ('norm' x) ('norm' v)
--     tol = 1e-3
-- @
--
-- <<docs/differential_equations/geodesic_hill_and_valley.svg>>
geodesicEquation
    :: (Double -> Vec2 -> Double) -- ^ Surface function \(f(t, \mathbf v)\)
    -> Double                     -- ^ Time \(t\)
    -> (Vec2, Vec2)               -- ^ \((\mathbf v, \dot{\mathbf v})\)
    -> (Vec2, Vec2)               -- ^ \((\dot{\mathbf v}, \ddot{\mathbf v})\)
geodesicEquation :: (Double -> Vec2 -> Double)
-> Double -> (Vec2, Vec2) -> (Vec2, Vec2)
geodesicEquation Double -> Vec2 -> Double
f Double
t (Vec2
v, v' :: Vec2
v'@(Vec2 Double
x' Double
y')) =
    ( Vec2
v'
    , Double -> Double -> Vec2
Vec2
        (-Dim -> Dim -> Double
c'x__V Dim
X Dim
XDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
x'Double -> Integer -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Dim -> Dim -> Double
c'x__V Dim
X Dim
YDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
x'Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
y' Double -> Double -> Double
forall a. Num a => a -> a -> a
-Dim -> Dim -> Double
c'x__V Dim
Y Dim
YDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
y'Double -> Integer -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2)
        (-Dim -> Dim -> Double
c'y__V Dim
X Dim
XDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
x'Double -> Integer -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Dim -> Dim -> Double
c'y__V Dim
X Dim
YDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
x'Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
y' Double -> Double -> Double
forall a. Num a => a -> a -> a
-Dim -> Dim -> Double
c'y__V Dim
Y Dim
YDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
y'Double -> Integer -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2)
    )
  where
    h :: Double
h = Double
1e-3

    -- Offsets for first derivatives at our current position
    vXH :: Vec2
vXH = Vec2
v Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
+. Double -> Double -> Vec2
Vec2 Double
h Double
0
    vYH :: Vec2
vYH = Vec2
v Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
+. Double -> Double -> Vec2
Vec2 Double
0 Double
h

    -- Offsets for second derivatives at our current position
    vXHXH :: Vec2
vXHXH = Vec2
vXH Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
+. Double -> Double -> Vec2
Vec2 Double
h Double
0
    vXHYH :: Vec2
vXHYH = Vec2
vYH Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
+. Double -> Double -> Vec2
Vec2 Double
h Double
0 -- = vYHXH – Vector addition commutativity saves us another call to f!
    vYHYH :: Vec2
vYHYH = Vec2
vYH Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
+. Double -> Double -> Vec2
Vec2 Double
0 Double
h

    -- Function application sharing
    ftv :: Double
ftv = Double -> Vec2 -> Double
f Double
t Vec2
v
    ftvXH :: Double
ftvXH = Double -> Vec2 -> Double
f Double
t Vec2
vXH
    ftvYH :: Double
ftvYH = Double -> Vec2 -> Double
f Double
t Vec2
vYH
    ftvXHYH :: Double
ftvXHYH = Double -> Vec2 -> Double
f Double
t Vec2
vXHYH

    -- First derivatives, applied to the offsets, for the second derivatives
    fdxvXH :: Double
fdxvXH = (Double -> Vec2 -> Double
f Double
t Vec2
vXHXH Double -> Double -> Double
forall v. VectorSpace v => v -> v -> v
-. Double
ftvXH) Double -> Double -> Double
forall v. VectorSpace v => v -> Double -> v
/. Double
h
    fdxvYH :: Double
fdxvYH = (Double
ftvXHYH   Double -> Double -> Double
forall v. VectorSpace v => v -> v -> v
-. Double
ftvYH) Double -> Double -> Double
forall v. VectorSpace v => v -> Double -> v
/. Double
h
    fdyvXH :: Double
fdyvXH = (Double
ftvXHYH   Double -> Double -> Double
forall v. VectorSpace v => v -> v -> v
-. Double
ftvXH) Double -> Double -> Double
forall v. VectorSpace v => v -> Double -> v
/. Double
h
    fdyvYH :: Double
fdyvYH = (Double -> Vec2 -> Double
f Double
t Vec2
vYHYH Double -> Double -> Double
forall v. VectorSpace v => v -> v -> v
-. Double
ftvYH) Double -> Double -> Double
forall v. VectorSpace v => v -> Double -> v
/. Double
h

    -- First derivatives applied at our current position
    fdxV :: Double
fdxV = (Double
ftvXH Double -> Double -> Double
forall v. VectorSpace v => v -> v -> v
-. Double
ftv) Double -> Double -> Double
forall v. VectorSpace v => v -> Double -> v
/. Double
h
    fdyV :: Double
fdyV = (Double
ftvYH Double -> Double -> Double
forall v. VectorSpace v => v -> v -> v
-. Double
ftv) Double -> Double -> Double
forall v. VectorSpace v => v -> Double -> v
/. Double
h

    -- Inverse metric g^{ab}
    (Double
g'x'x, Double
g'x'y, Double
g'y'x, Double
g'y'y) =
        let denominator :: Double
denominator = (Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
fdxVDouble -> Integer -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
fdyVDouble -> Integer -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2)
        in ( (Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
fdyVDouble -> Integer -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2)   Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
denominator
           , -(Double
fdxVDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
fdyV) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
denominator
           , Double
g'x'y
           , (Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
fdxVDouble -> Integer -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2)   Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
denominator
        )

    -- Derivative of the metric g_{ab,c}
    g__d_V :: Dim -> Dim -> Dim -> Double
g__d_V Dim
X Dim
X Dim
X = Double
g_x_xd_xV
    g__d_V Dim
X Dim
X Dim
Y = Double
g_x_xd_yV
    g__d_V Dim
X Dim
Y Dim
X = Double
g_x_yd_xV
    g__d_V Dim
X Dim
Y Dim
Y = Double
g_x_yd_yV
    g__d_V Dim
Y Dim
X Dim
X = Double
g_y_xd_xV
    g__d_V Dim
Y Dim
X Dim
Y = Double
g_y_xd_yV
    g__d_V Dim
Y Dim
Y Dim
X = Double
g_y_yd_xV
    g__d_V Dim
Y Dim
Y Dim
Y = Double
g_y_yd_yV

    -- Derivative of the metric g_{ab,c} at our current position
    g_x_xd_xV :: Double
g_x_xd_xV = ((Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
fdxvXHDouble -> Integer -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2)    Double -> Double -> Double
forall v. VectorSpace v => v -> v -> v
-. (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
fdxVDouble -> Integer -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2))  Double -> Double -> Double
forall v. VectorSpace v => v -> Double -> v
/. Double
h
    g_x_xd_yV :: Double
g_x_xd_yV = ((Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
fdxvYHDouble -> Integer -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2)    Double -> Double -> Double
forall v. VectorSpace v => v -> v -> v
-. (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
fdxVDouble -> Integer -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2))  Double -> Double -> Double
forall v. VectorSpace v => v -> Double -> v
/. Double
h
    g_x_yd_xV :: Double
g_x_yd_xV = ((Double
fdxvXH Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
fdyvXH) Double -> Double -> Double
forall v. VectorSpace v => v -> v -> v
-. (Double
fdxV Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
fdyV)) Double -> Double -> Double
forall v. VectorSpace v => v -> Double -> v
/. Double
h
    g_x_yd_yV :: Double
g_x_yd_yV = ((Double
fdxvYH Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
fdyvYH) Double -> Double -> Double
forall v. VectorSpace v => v -> v -> v
-. (Double
fdxV Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
fdyV)) Double -> Double -> Double
forall v. VectorSpace v => v -> Double -> v
/. Double
h
    g_y_xd_xV :: Double
g_y_xd_xV = ((Double
fdxvXH Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
fdyvXH) Double -> Double -> Double
forall v. VectorSpace v => v -> v -> v
-. (Double
fdxV Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
fdyV)) Double -> Double -> Double
forall v. VectorSpace v => v -> Double -> v
/. Double
h
    g_y_xd_yV :: Double
g_y_xd_yV = ((Double
fdxvYH Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
fdyvYH) Double -> Double -> Double
forall v. VectorSpace v => v -> v -> v
-. (Double
fdxV Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
fdyV)) Double -> Double -> Double
forall v. VectorSpace v => v -> Double -> v
/. Double
h
    g_y_yd_xV :: Double
g_y_yd_xV = ((Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
fdyvXHDouble -> Integer -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2)    Double -> Double -> Double
forall v. VectorSpace v => v -> v -> v
-. (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
fdyVDouble -> Integer -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2))  Double -> Double -> Double
forall v. VectorSpace v => v -> Double -> v
/. Double
h
    g_y_yd_yV :: Double
g_y_yd_yV = ((Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
fdyvYHDouble -> Integer -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2)    Double -> Double -> Double
forall v. VectorSpace v => v -> v -> v
-. (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
fdyVDouble -> Integer -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2))  Double -> Double -> Double
forall v. VectorSpace v => v -> Double -> v
/. Double
h

    -- Christoffel symbols, \Gamma^i_{kl} = \frac12 g^{im} (g_{mk,l}+g_{ml,k}-g_{kl,m})
    c'x__V :: Dim -> Dim -> Double
c'x__V Dim
k Dim
l = Double
0.5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
g'x'x Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Dim -> Dim -> Dim -> Double
g__d_V Dim
X Dim
k Dim
l Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Dim -> Dim -> Dim -> Double
g__d_V Dim
X Dim
l Dim
k Double -> Double -> Double
forall a. Num a => a -> a -> a
- Dim -> Dim -> Dim -> Double
g__d_V Dim
k Dim
l Dim
X) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
g'x'y Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Dim -> Dim -> Dim -> Double
g__d_V Dim
Y Dim
k Dim
l Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Dim -> Dim -> Dim -> Double
g__d_V Dim
Y Dim
l Dim
k Double -> Double -> Double
forall a. Num a => a -> a -> a
- Dim -> Dim -> Dim -> Double
g__d_V Dim
k Dim
l Dim
Y))
    c'y__V :: Dim -> Dim -> Double
c'y__V Dim
k Dim
l = Double
0.5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
g'y'x Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Dim -> Dim -> Dim -> Double
g__d_V Dim
X Dim
k Dim
l Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Dim -> Dim -> Dim -> Double
g__d_V Dim
X Dim
l Dim
k Double -> Double -> Double
forall a. Num a => a -> a -> a
- Dim -> Dim -> Dim -> Double
g__d_V Dim
k Dim
l Dim
X) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
g'y'y Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Dim -> Dim -> Dim -> Double
g__d_V Dim
Y Dim
k Dim
l Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Dim -> Dim -> Dim -> Double
g__d_V Dim
Y Dim
l Dim
k Double -> Double -> Double
forall a. Num a => a -> a -> a
- Dim -> Dim -> Dim -> Double
g__d_V Dim
k Dim
l Dim
Y))

data Dim = X | Y