{-# LANGUAGE RebindableSyntax #-}

-- | Riemann curvature from a Levi-Civita connection.
--
-- Following the classical identity (Albert arXiv:2312.02664 Def. 6)
--
-- > R^ρ_{σμν} = ∂_μ Γ^ρ_{νσ} − ∂_ν Γ^ρ_{μσ} + Γ^ρ_{μλ} Γ^λ_{νσ} − Γ^ρ_{νλ} Γ^λ_{μσ}
--
-- built on 'christoffel2D' for 2-D coordinate metrics, plus a diagonal-metric
-- path for the 4-D Schwarzschild vacuum oracle.
module Circuit.Diff.Curvature
  ( -- * 2-D Riemann / Ricci
    Gamma2D,
    gamma2DAt,
    riemann2D,
    ricci2D,
    ricciScalar2D,

    -- * 2-D metrics
    sphereMetricLower,
    sphereMetricRaise,

    -- * Diagonal n-D (Schwarzschild)
    DiagonalMetric (..),
    schwarzschildMetric,
    gammaDiagonal,
    riemannDiagonal,
    ricciDiagonal,
  )
where

import Circuit.Diff (Diff (..), Diff', runDiff)
import Circuit.Diff.Chart
  ( christoffel2D,
    raise2D,
  )
import NumHask.Prelude

-- | Packed Christoffel symbols @Γ^c_{ab}@ for a 2-D chart (c,a,b ∈ {0,1}).
type Gamma2D a = (a, a, a, a, a, a, a, a)

-- | @Γ^c_{ab}@ lookup.
gammaComp :: Gamma2D a -> Int -> Int -> Int -> a
gammaComp :: forall a. Gamma2D a -> Int -> Int -> Int -> a
gammaComp (a
g000, a
g001, a
g010, a
g011, a
g100, a
g101, a
g110, a
g111) Int
c Int
a Int
b =
  case (Int
c, Int
a, Int
b) of
    (Int
0, Int
0, Int
0) -> a
g000
    (Int
0, Int
0, Int
1) -> a
g001
    (Int
0, Int
1, Int
0) -> a
g010
    (Int
0, Int
1, Int
1) -> a
g011
    (Int
1, Int
0, Int
0) -> a
g100
    (Int
1, Int
0, Int
1) -> a
g101
    (Int
1, Int
1, Int
0) -> a
g110
    (Int
1, Int
1, Int
1) -> a
g111
    (Int, Int, Int)
_ -> [Char] -> a
forall a. HasCallStack => [Char] -> a
error [Char]
"gammaComp: index out of range"

gamma2DAt ::
  (Field a) =>
  Diff' ((a, a), (a, a)) (a, a) ->
  Diff' ((a, a), (a, a)) (a, a) ->
  (a, a) ->
  Gamma2D a
gamma2DAt :: forall a.
Field a =>
Diff' ((a, a), (a, a)) (a, a)
-> Diff' ((a, a), (a, a)) (a, a) -> (a, a) -> Gamma2D a
gamma2DAt = Diff' ((a, a), (a, a)) (a, a)
-> Diff' ((a, a), (a, a)) (a, a)
-> (a, a)
-> (a, a, a, a, a, a, a, a)
forall a.
Field a =>
Diff' ((a, a), (a, a)) (a, a)
-> Diff' ((a, a), (a, a)) (a, a) -> (a, a) -> Gamma2D a
christoffel2D

-- | Central difference step for ∂Γ.
fdStep :: (Field a, FromInteger a) => a
fdStep :: forall a. (Field a, FromInteger a) => a
fdStep = a
forall a. Multiplicative a => a
one a -> a -> a
forall a. Divisive a => a -> a -> a
/ Integer -> a
forall a. FromInteger a => Integer -> a
fromInteger (Integer
10000000 :: Integer)

-- | Basis step along coordinate @i@.
bump2 :: (Additive a) => Int -> a -> (a, a) -> (a, a)
bump2 :: forall a. Additive a => Int -> a -> (a, a) -> (a, a)
bump2 Int
0 a
h (a
x0, a
x1) = (a
x0 a -> a -> a
forall a. Additive a => a -> a -> a
+ a
h, a
x1)
bump2 Int
1 a
h (a
x0, a
x1) = (a
x0, a
x1 a -> a -> a
forall a. Additive a => a -> a -> a
+ a
h)
bump2 Int
_ a
_ (a, a)
_ = [Char] -> (a, a)
forall a. HasCallStack => [Char] -> a
error [Char]
"bump2: bad index"

-- | @∂_μ Γ^ρ_{ab}@ by central differences of 'christoffel2D'.
partialGamma2D ::
  (Field a, FromInteger a) =>
  Diff' ((a, a), (a, a)) (a, a) ->
  Diff' ((a, a), (a, a)) (a, a) ->
  (a, a) ->
  Int ->
  Int ->
  Int ->
  Int ->
  a
partialGamma2D :: forall a.
(Field a, FromInteger a) =>
Diff' ((a, a), (a, a)) (a, a)
-> Diff' ((a, a), (a, a)) (a, a)
-> (a, a)
-> Int
-> Int
-> Int
-> Int
-> a
partialGamma2D Diff' ((a, a), (a, a)) (a, a)
lower Diff' ((a, a), (a, a)) (a, a)
raise (a, a)
x Int
mu Int
rho Int
a Int
b =
  let h :: a
h = a
forall a. (Field a, FromInteger a) => a
fdStep
      gPlus :: Gamma2D a
gPlus = Diff' ((a, a), (a, a)) (a, a)
-> Diff' ((a, a), (a, a)) (a, a) -> (a, a) -> Gamma2D a
forall a.
Field a =>
Diff' ((a, a), (a, a)) (a, a)
-> Diff' ((a, a), (a, a)) (a, a) -> (a, a) -> Gamma2D a
gamma2DAt Diff' ((a, a), (a, a)) (a, a)
lower Diff' ((a, a), (a, a)) (a, a)
raise (Int -> a -> (a, a) -> (a, a)
forall a. Additive a => Int -> a -> (a, a) -> (a, a)
bump2 Int
mu a
h (a, a)
x)
      gMinus :: Gamma2D a
gMinus = Diff' ((a, a), (a, a)) (a, a)
-> Diff' ((a, a), (a, a)) (a, a) -> (a, a) -> Gamma2D a
forall a.
Field a =>
Diff' ((a, a), (a, a)) (a, a)
-> Diff' ((a, a), (a, a)) (a, a) -> (a, a) -> Gamma2D a
gamma2DAt Diff' ((a, a), (a, a)) (a, a)
lower Diff' ((a, a), (a, a)) (a, a)
raise (Int -> a -> (a, a) -> (a, a)
forall a. Additive a => Int -> a -> (a, a) -> (a, a)
bump2 Int
mu (a -> a
forall a. Subtractive a => a -> a
negate a
h) (a, a)
x)
      twoH :: a
twoH = a
h a -> a -> a
forall a. Additive a => a -> a -> a
+ a
h
   in (Gamma2D a -> Int -> Int -> Int -> a
forall a. Gamma2D a -> Int -> Int -> Int -> a
gammaComp Gamma2D a
gPlus Int
rho Int
a Int
b a -> a -> a
forall a. Subtractive a => a -> a -> a
- Gamma2D a -> Int -> Int -> Int -> a
forall a. Gamma2D a -> Int -> Int -> Int -> a
gammaComp Gamma2D a
gMinus Int
rho Int
a Int
b) a -> a -> a
forall a. Divisive a => a -> a -> a
/ a
twoH

-- | Riemann component @R^ρ_{σμν}@ at a point (2-D).
--
-- > R^ρ_{σμν} = ∂_μ Γ^ρ_{νσ} − ∂_ν Γ^ρ_{μσ}
-- >           + Γ^ρ_{μλ} Γ^λ_{νσ} − Γ^ρ_{νλ} Γ^λ_{μσ}
riemann2D ::
  (Field a, FromInteger a) =>
  Diff' ((a, a), (a, a)) (a, a) ->
  Diff' ((a, a), (a, a)) (a, a) ->
  (a, a) ->
  Int ->
  Int ->
  Int ->
  Int ->
  a
riemann2D :: forall a.
(Field a, FromInteger a) =>
Diff' ((a, a), (a, a)) (a, a)
-> Diff' ((a, a), (a, a)) (a, a)
-> (a, a)
-> Int
-> Int
-> Int
-> Int
-> a
riemann2D Diff' ((a, a), (a, a)) (a, a)
lower Diff' ((a, a), (a, a)) (a, a)
raise (a, a)
x Int
rho Int
sigma Int
mu Int
nu =
  let g :: Gamma2D a
g = Diff' ((a, a), (a, a)) (a, a)
-> Diff' ((a, a), (a, a)) (a, a) -> (a, a) -> Gamma2D a
forall a.
Field a =>
Diff' ((a, a), (a, a)) (a, a)
-> Diff' ((a, a), (a, a)) (a, a) -> (a, a) -> Gamma2D a
gamma2DAt Diff' ((a, a), (a, a)) (a, a)
lower Diff' ((a, a), (a, a)) (a, a)
raise (a, a)
x
      dMu :: a
dMu = Diff' ((a, a), (a, a)) (a, a)
-> Diff' ((a, a), (a, a)) (a, a)
-> (a, a)
-> Int
-> Int
-> Int
-> Int
-> a
forall a.
(Field a, FromInteger a) =>
Diff' ((a, a), (a, a)) (a, a)
-> Diff' ((a, a), (a, a)) (a, a)
-> (a, a)
-> Int
-> Int
-> Int
-> Int
-> a
partialGamma2D Diff' ((a, a), (a, a)) (a, a)
lower Diff' ((a, a), (a, a)) (a, a)
raise (a, a)
x Int
mu Int
rho Int
nu Int
sigma
      dNu :: a
dNu = Diff' ((a, a), (a, a)) (a, a)
-> Diff' ((a, a), (a, a)) (a, a)
-> (a, a)
-> Int
-> Int
-> Int
-> Int
-> a
forall a.
(Field a, FromInteger a) =>
Diff' ((a, a), (a, a)) (a, a)
-> Diff' ((a, a), (a, a)) (a, a)
-> (a, a)
-> Int
-> Int
-> Int
-> Int
-> a
partialGamma2D Diff' ((a, a), (a, a)) (a, a)
lower Diff' ((a, a), (a, a)) (a, a)
raise (a, a)
x Int
nu Int
rho Int
mu Int
sigma
      quad :: a
quad =
        [a] -> a
forall a (f :: * -> *). (Additive a, Foldable f) => f a -> a
sum
          [ Gamma2D a -> Int -> Int -> Int -> a
forall a. Gamma2D a -> Int -> Int -> Int -> a
gammaComp Gamma2D a
g Int
rho Int
mu Int
lam
              a -> a -> a
forall a. Multiplicative a => a -> a -> a
* Gamma2D a -> Int -> Int -> Int -> a
forall a. Gamma2D a -> Int -> Int -> Int -> a
gammaComp Gamma2D a
g Int
lam Int
nu Int
sigma
              a -> a -> a
forall a. Subtractive a => a -> a -> a
- Gamma2D a -> Int -> Int -> Int -> a
forall a. Gamma2D a -> Int -> Int -> Int -> a
gammaComp Gamma2D a
g Int
rho Int
nu Int
lam
              a -> a -> a
forall a. Multiplicative a => a -> a -> a
* Gamma2D a -> Int -> Int -> Int -> a
forall a. Gamma2D a -> Int -> Int -> Int -> a
gammaComp Gamma2D a
g Int
lam Int
mu Int
sigma
          | Int
lam <- [Int
0, Int
1]
          ]
   in a
dMu a -> a -> a
forall a. Subtractive a => a -> a -> a
- a
dNu a -> a -> a
forall a. Additive a => a -> a -> a
+ a
quad

-- | Ricci @R_{σν} = R^ρ_{σρν}@ (sum on ρ).
ricci2D ::
  (Field a, FromInteger a) =>
  Diff' ((a, a), (a, a)) (a, a) ->
  Diff' ((a, a), (a, a)) (a, a) ->
  (a, a) ->
  Int ->
  Int ->
  a
ricci2D :: forall a.
(Field a, FromInteger a) =>
Diff' ((a, a), (a, a)) (a, a)
-> Diff' ((a, a), (a, a)) (a, a) -> (a, a) -> Int -> Int -> a
ricci2D Diff' ((a, a), (a, a)) (a, a)
lower Diff' ((a, a), (a, a)) (a, a)
raise (a, a)
x Int
sigma Int
nu =
  [a] -> a
forall a (f :: * -> *). (Additive a, Foldable f) => f a -> a
sum [Diff' ((a, a), (a, a)) (a, a)
-> Diff' ((a, a), (a, a)) (a, a)
-> (a, a)
-> Int
-> Int
-> Int
-> Int
-> a
forall a.
(Field a, FromInteger a) =>
Diff' ((a, a), (a, a)) (a, a)
-> Diff' ((a, a), (a, a)) (a, a)
-> (a, a)
-> Int
-> Int
-> Int
-> Int
-> a
riemann2D Diff' ((a, a), (a, a)) (a, a)
lower Diff' ((a, a), (a, a)) (a, a)
raise (a, a)
x Int
rho Int
sigma Int
rho Int
nu | Int
rho <- [Int
0, Int
1]]

-- | Ricci scalar @R = g^{σν} R_{σν}@.
ricciScalar2D ::
  (Field a, FromInteger a) =>
  Diff' ((a, a), (a, a)) (a, a) ->
  Diff' ((a, a), (a, a)) (a, a) ->
  (a, a) ->
  a
ricciScalar2D :: forall a.
(Field a, FromInteger a) =>
Diff' ((a, a), (a, a)) (a, a)
-> Diff' ((a, a), (a, a)) (a, a) -> (a, a) -> a
ricciScalar2D Diff' ((a, a), (a, a)) (a, a)
lower Diff' ((a, a), (a, a)) (a, a)
raise (a, a)
x =
  let ((a, a)
gInv0, (a, a) -> ((a, a), (a, a))
_) = Diff' ((a, a), (a, a)) (a, a)
-> ((a, a), (a, a)) -> ((a, a), (a, a) -> ((a, a), (a, a)))
forall k (p :: k) a b. Diff p a b -> a -> (b, b -> a)
runDiff Diff' ((a, a), (a, a)) (a, a)
raise ((a, a)
x, (a
forall a. Multiplicative a => a
one, a
forall a. Additive a => a
zero))
      ((a, a)
gInv1, (a, a) -> ((a, a), (a, a))
_) = Diff' ((a, a), (a, a)) (a, a)
-> ((a, a), (a, a)) -> ((a, a), (a, a) -> ((a, a), (a, a)))
forall k (p :: k) a b. Diff p a b -> a -> (b, b -> a)
runDiff Diff' ((a, a), (a, a)) (a, a)
raise ((a, a)
x, (a
forall a. Additive a => a
zero, a
forall a. Multiplicative a => a
one))
      g00 :: a
g00 = (a, a) -> a
forall a b. (a, b) -> a
fst (a, a)
gInv0
      g01 :: a
g01 = (a, a) -> a
forall a b. (a, b) -> a
fst (a, a)
gInv1
      g10 :: a
g10 = (a, a) -> a
forall a b. (a, b) -> b
snd (a, a)
gInv0
      g11 :: a
g11 = (a, a) -> a
forall a b. (a, b) -> b
snd (a, a)
gInv1
      r00 :: a
r00 = Diff' ((a, a), (a, a)) (a, a)
-> Diff' ((a, a), (a, a)) (a, a) -> (a, a) -> Int -> Int -> a
forall a.
(Field a, FromInteger a) =>
Diff' ((a, a), (a, a)) (a, a)
-> Diff' ((a, a), (a, a)) (a, a) -> (a, a) -> Int -> Int -> a
ricci2D Diff' ((a, a), (a, a)) (a, a)
lower Diff' ((a, a), (a, a)) (a, a)
raise (a, a)
x Int
0 Int
0
      r01 :: a
r01 = Diff' ((a, a), (a, a)) (a, a)
-> Diff' ((a, a), (a, a)) (a, a) -> (a, a) -> Int -> Int -> a
forall a.
(Field a, FromInteger a) =>
Diff' ((a, a), (a, a)) (a, a)
-> Diff' ((a, a), (a, a)) (a, a) -> (a, a) -> Int -> Int -> a
ricci2D Diff' ((a, a), (a, a)) (a, a)
lower Diff' ((a, a), (a, a)) (a, a)
raise (a, a)
x Int
0 Int
1
      r10 :: a
r10 = Diff' ((a, a), (a, a)) (a, a)
-> Diff' ((a, a), (a, a)) (a, a) -> (a, a) -> Int -> Int -> a
forall a.
(Field a, FromInteger a) =>
Diff' ((a, a), (a, a)) (a, a)
-> Diff' ((a, a), (a, a)) (a, a) -> (a, a) -> Int -> Int -> a
ricci2D Diff' ((a, a), (a, a)) (a, a)
lower Diff' ((a, a), (a, a)) (a, a)
raise (a, a)
x Int
1 Int
0
      r11 :: a
r11 = Diff' ((a, a), (a, a)) (a, a)
-> Diff' ((a, a), (a, a)) (a, a) -> (a, a) -> Int -> Int -> a
forall a.
(Field a, FromInteger a) =>
Diff' ((a, a), (a, a)) (a, a)
-> Diff' ((a, a), (a, a)) (a, a) -> (a, a) -> Int -> Int -> a
ricci2D Diff' ((a, a), (a, a)) (a, a)
lower Diff' ((a, a), (a, a)) (a, a)
raise (a, a)
x Int
1 Int
1
   in a
g00 a -> a -> a
forall a. Multiplicative a => a -> a -> a
* a
r00 a -> a -> a
forall a. Additive a => a -> a -> a
+ a
g01 a -> a -> a
forall a. Multiplicative a => a -> a -> a
* a
r01 a -> a -> a
forall a. Additive a => a -> a -> a
+ a
g10 a -> a -> a
forall a. Multiplicative a => a -> a -> a
* a
r10 a -> a -> a
forall a. Additive a => a -> a -> a
+ a
g11 a -> a -> a
forall a. Multiplicative a => a -> a -> a
* a
r11

--------------------------------------------------------------------------------
-- Unit sphere metric (θ, φ):  g = diag(1, sin²θ)
--------------------------------------------------------------------------------

-- | Unit 2-sphere metric in spherical coordinates @(θ, φ)@.
--
-- @ds² = dθ² + sin²θ dφ²@.  Analytic Ricci scalar is @2@.
sphereMetricLower :: (TrigField a) => Diff' ((a, a), (a, a)) (a, a)
sphereMetricLower :: forall a. TrigField a => Diff' ((a, a), (a, a)) (a, a)
sphereMetricLower = (((a, a), (a, a)) -> ((a, a), (a, a) -> ((a, a), (a, a))))
-> Diff () ((a, a), (a, a)) (a, a)
forall k (p :: k) a b. (a -> (b, b -> a)) -> Diff p a b
Diff ((((a, a), (a, a)) -> ((a, a), (a, a) -> ((a, a), (a, a))))
 -> Diff () ((a, a), (a, a)) (a, a))
-> (((a, a), (a, a)) -> ((a, a), (a, a) -> ((a, a), (a, a))))
-> Diff () ((a, a), (a, a)) (a, a)
forall a b. (a -> b) -> a -> b
$ \((a
theta, a
_), (a
vth, a
vph)) ->
  let s :: a
s = a -> a
forall a. TrigField a => a -> a
sin a
theta
      ss :: a
ss = a
s a -> a -> a
forall a. Multiplicative a => a -> a -> a
* a
s
      tw :: a
tw = a
forall a. Multiplicative a => a
one a -> a -> a
forall a. Additive a => a -> a -> a
+ a
forall a. Multiplicative a => a
one
   in ( (a
vth, a
ss a -> a -> a
forall a. Multiplicative a => a -> a -> a
* a
vph),
        \(a
dcth, a
dcph) ->
          ( (a
tw a -> a -> a
forall a. Multiplicative a => a -> a -> a
* a
s a -> a -> a
forall a. Multiplicative a => a -> a -> a
* a -> a
forall a. TrigField a => a -> a
cos a
theta a -> a -> a
forall a. Multiplicative a => a -> a -> a
* a
vph a -> a -> a
forall a. Multiplicative a => a -> a -> a
* a
dcph, a
forall a. Additive a => a
zero),
            (a
dcth, a
ss a -> a -> a
forall a. Multiplicative a => a -> a -> a
* a
dcph)
          )
      )

sphereMetricRaise :: (TrigField a) => Diff' ((a, a), (a, a)) (a, a)
sphereMetricRaise :: forall a. TrigField a => Diff' ((a, a), (a, a)) (a, a)
sphereMetricRaise = Diff' ((a, a), (a, a)) (a, a) -> Diff' ((a, a), (a, a)) (a, a)
forall a.
Field a =>
Diff' ((a, a), (a, a)) (a, a) -> Diff' ((a, a), (a, a)) (a, a)
raise2D Diff' ((a, a), (a, a)) (a, a)
forall a. TrigField a => Diff' ((a, a), (a, a)) (a, a)
sphereMetricLower

--------------------------------------------------------------------------------
-- Diagonal metrics (Schwarzschild)
--------------------------------------------------------------------------------

-- | A diagonal metric given by its diagonal components @g_ii(x)@ and their
-- first partials @∂_j g_ii(x)@ (no sum).
data DiagonalMetric a = DiagonalMetric
  { forall a. DiagonalMetric a -> Int
dmDim :: Int,
    forall a. DiagonalMetric a -> [a] -> Int -> a
dmG :: [a] -> Int -> a,
    forall a. DiagonalMetric a -> [a] -> Int -> Int -> a
dmDG :: [a] -> Int -> Int -> a
  }

-- | Schwarzschild metric outside the horizon (@r > 2M@).
--
-- Coordinates @(t, r, θ, φ)@; @rs = 2M@ is the Schwarzschild radius.
--
-- > g = diag( −(1−rs/r), 1/(1−rs/r), r², r² sin²θ )
schwarzschildMetric :: (TrigField a) => a -> DiagonalMetric a
schwarzschildMetric :: forall a. TrigField a => a -> DiagonalMetric a
schwarzschildMetric a
rs =
  DiagonalMetric
    { dmDim :: Int
dmDim = Int
4,
      dmG :: [a] -> Int -> a
dmG = \[a]
xs Int
i ->
        let r :: a
r = [a]
xs [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
1
            th :: a
th = [a]
xs [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
2
            f :: a
f = a
forall a. Multiplicative a => a
one a -> a -> a
forall a. Subtractive a => a -> a -> a
- a
rs a -> a -> a
forall a. Divisive a => a -> a -> a
/ a
r
         in case Int
i of
              Int
0 -> a -> a
forall a. Subtractive a => a -> a
negate a
f
              Int
1 -> a -> a
forall a. Divisive a => a -> a
recip a
f
              Int
2 -> a
r a -> a -> a
forall a. Multiplicative a => a -> a -> a
* a
r
              Int
3 -> a
r a -> a -> a
forall a. Multiplicative a => a -> a -> a
* a
r a -> a -> a
forall a. Multiplicative a => a -> a -> a
* a -> a
forall a. TrigField a => a -> a
sin a
th a -> a -> a
forall a. Multiplicative a => a -> a -> a
* a -> a
forall a. TrigField a => a -> a
sin a
th
              Int
_ -> [Char] -> a
forall a. HasCallStack => [Char] -> a
error [Char]
"schwarzschildMetric: bad index",
      dmDG :: [a] -> Int -> Int -> a
dmDG = \[a]
xs Int
j Int
i ->
        let r :: a
r = [a]
xs [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
1
            th :: a
th = [a]
xs [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
2
            f :: a
f = a
forall a. Multiplicative a => a
one a -> a -> a
forall a. Subtractive a => a -> a -> a
- a
rs a -> a -> a
forall a. Divisive a => a -> a -> a
/ a
r
            -- ∂_r f = rs/r²
            dFdr :: a
dFdr = a
rs a -> a -> a
forall a. Divisive a => a -> a -> a
/ (a
r a -> a -> a
forall a. Multiplicative a => a -> a -> a
* a
r)
         in case (Int
j, Int
i) of
              -- only r-derivatives of g_tt, g_rr; θ-derivative of g_φφ
              (Int
1, Int
0) -> a -> a
forall a. Subtractive a => a -> a
negate a
dFdr
              (Int
1, Int
1) -> a -> a
forall a. Subtractive a => a -> a
negate a
dFdr a -> a -> a
forall a. Divisive a => a -> a -> a
/ (a
f a -> a -> a
forall a. Multiplicative a => a -> a -> a
* a
f)
              (Int
1, Int
2) -> (a
forall a. Multiplicative a => a
one a -> a -> a
forall a. Additive a => a -> a -> a
+ a
forall a. Multiplicative a => a
one) a -> a -> a
forall a. Multiplicative a => a -> a -> a
* a
r
              (Int
1, Int
3) -> (a
forall a. Multiplicative a => a
one a -> a -> a
forall a. Additive a => a -> a -> a
+ a
forall a. Multiplicative a => a
one) a -> a -> a
forall a. Multiplicative a => a -> a -> a
* a
r a -> a -> a
forall a. Multiplicative a => a -> a -> a
* a -> a
forall a. TrigField a => a -> a
sin a
th a -> a -> a
forall a. Multiplicative a => a -> a -> a
* a -> a
forall a. TrigField a => a -> a
sin a
th
              (Int
2, Int
3) -> (a
forall a. Multiplicative a => a
one a -> a -> a
forall a. Additive a => a -> a -> a
+ a
forall a. Multiplicative a => a
one) a -> a -> a
forall a. Multiplicative a => a -> a -> a
* a
r a -> a -> a
forall a. Multiplicative a => a -> a -> a
* a
r a -> a -> a
forall a. Multiplicative a => a -> a -> a
* a -> a
forall a. TrigField a => a -> a
sin a
th a -> a -> a
forall a. Multiplicative a => a -> a -> a
* a -> a
forall a. TrigField a => a -> a
cos a
th
              (Int, Int)
_ -> a
forall a. Additive a => a
zero
    }

-- | Inverse diagonal: @g^{ii} = 1/g_ii@.
gInvDiag :: (Field a) => DiagonalMetric a -> [a] -> Int -> a
gInvDiag :: forall a. Field a => DiagonalMetric a -> [a] -> Int -> a
gInvDiag DiagonalMetric a
m [a]
x Int
i = a -> a
forall a. Divisive a => a -> a
recip (DiagonalMetric a -> [a] -> Int -> a
forall a. DiagonalMetric a -> [a] -> Int -> a
dmG DiagonalMetric a
m [a]
x Int
i)

-- | Christoffel @Γ^c_{ab}@ for a diagonal metric.
--
-- > Γ^c_{ab} = ½ g^{cc} (∂_a g_{bc} + ∂_b g_{ac} − ∂_c g_{ab})  (no sum on c)
gammaDiagonal :: (Field a) => DiagonalMetric a -> [a] -> Int -> Int -> Int -> a
gammaDiagonal :: forall a.
Field a =>
DiagonalMetric a -> [a] -> Int -> Int -> Int -> a
gammaDiagonal DiagonalMetric a
m [a]
x Int
c Int
a Int
b =
  let gInv :: a
gInv = DiagonalMetric a -> [a] -> Int -> a
forall a. Field a => DiagonalMetric a -> [a] -> Int -> a
gInvDiag DiagonalMetric a
m [a]
x Int
c
      -- g_ab is zero unless a==b for diagonal metrics
      dg :: Int -> Int -> Int -> a
dg Int
a' Int
b' Int
c' =
        if Int
a' Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
b'
          then DiagonalMetric a -> [a] -> Int -> Int -> a
forall a. DiagonalMetric a -> [a] -> Int -> Int -> a
dmDG DiagonalMetric a
m [a]
x Int
c' Int
a'
          else a
forall a. Additive a => a
zero
      term :: a
term = Int -> Int -> Int -> a
dg Int
b Int
c Int
a a -> a -> a
forall a. Additive a => a -> a -> a
+ Int -> Int -> Int -> a
dg Int
a Int
c Int
b a -> a -> a
forall a. Subtractive a => a -> a -> a
- Int -> Int -> Int -> a
dg Int
a Int
b Int
c
   in a
forall a. (Additive a, Divisive a) => a
half a -> a -> a
forall a. Multiplicative a => a -> a -> a
* a
gInv a -> a -> a
forall a. Multiplicative a => a -> a -> a
* a
term

-- | @∂_μ Γ^c_{ab}@ by central differences on coordinate @μ@.
partialGammaDiagonal ::
  (Field a, FromInteger a) =>
  DiagonalMetric a ->
  [a] ->
  Int ->
  Int ->
  Int ->
  Int ->
  a
partialGammaDiagonal :: forall a.
(Field a, FromInteger a) =>
DiagonalMetric a -> [a] -> Int -> Int -> Int -> Int -> a
partialGammaDiagonal DiagonalMetric a
m [a]
x Int
mu Int
c Int
a Int
b =
  let h :: a
h = a
forall a. (Field a, FromInteger a) => a
fdStep
      n :: Int
n = DiagonalMetric a -> Int
forall a. DiagonalMetric a -> Int
dmDim DiagonalMetric a
m
      bump :: a -> [a]
bump a
j =
        [ if Int
k Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
mu then ([a]
x [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
k) a -> a -> a
forall a. Additive a => a -> a -> a
+ a
j a -> a -> a
forall a. Multiplicative a => a -> a -> a
* a
h else [a]
x [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
k
        | Int
k <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Subtractive a => a -> a -> a
- Int
1]
        ]
      twoH :: a
twoH = a
h a -> a -> a
forall a. Additive a => a -> a -> a
+ a
h
   in ( DiagonalMetric a -> [a] -> Int -> Int -> Int -> a
forall a.
Field a =>
DiagonalMetric a -> [a] -> Int -> Int -> Int -> a
gammaDiagonal DiagonalMetric a
m (a -> [a]
bump a
forall a. Multiplicative a => a
one) Int
c Int
a Int
b
          a -> a -> a
forall a. Subtractive a => a -> a -> a
- DiagonalMetric a -> [a] -> Int -> Int -> Int -> a
forall a.
Field a =>
DiagonalMetric a -> [a] -> Int -> Int -> Int -> a
gammaDiagonal DiagonalMetric a
m (a -> [a]
bump (a -> a
forall a. Subtractive a => a -> a
negate a
forall a. Multiplicative a => a
one)) Int
c Int
a Int
b
      )
        a -> a -> a
forall a. Divisive a => a -> a -> a
/ a
twoH

-- | Riemann @R^ρ_{σμν}@ for a diagonal metric.
riemannDiagonal ::
  (Field a, FromInteger a) =>
  DiagonalMetric a ->
  [a] ->
  Int ->
  Int ->
  Int ->
  Int ->
  a
riemannDiagonal :: forall a.
(Field a, FromInteger a) =>
DiagonalMetric a -> [a] -> Int -> Int -> Int -> Int -> a
riemannDiagonal DiagonalMetric a
m [a]
x Int
rho Int
sigma Int
mu Int
nu =
  let n :: Int
n = DiagonalMetric a -> Int
forall a. DiagonalMetric a -> Int
dmDim DiagonalMetric a
m
      dMu :: a
dMu = DiagonalMetric a -> [a] -> Int -> Int -> Int -> Int -> a
forall a.
(Field a, FromInteger a) =>
DiagonalMetric a -> [a] -> Int -> Int -> Int -> Int -> a
partialGammaDiagonal DiagonalMetric a
m [a]
x Int
mu Int
rho Int
nu Int
sigma
      dNu :: a
dNu = DiagonalMetric a -> [a] -> Int -> Int -> Int -> Int -> a
forall a.
(Field a, FromInteger a) =>
DiagonalMetric a -> [a] -> Int -> Int -> Int -> Int -> a
partialGammaDiagonal DiagonalMetric a
m [a]
x Int
nu Int
rho Int
mu Int
sigma
      quad :: a
quad =
        [a] -> a
forall a (f :: * -> *). (Additive a, Foldable f) => f a -> a
sum
          [ DiagonalMetric a -> [a] -> Int -> Int -> Int -> a
forall a.
Field a =>
DiagonalMetric a -> [a] -> Int -> Int -> Int -> a
gammaDiagonal DiagonalMetric a
m [a]
x Int
rho Int
mu Int
lam
              a -> a -> a
forall a. Multiplicative a => a -> a -> a
* DiagonalMetric a -> [a] -> Int -> Int -> Int -> a
forall a.
Field a =>
DiagonalMetric a -> [a] -> Int -> Int -> Int -> a
gammaDiagonal DiagonalMetric a
m [a]
x Int
lam Int
nu Int
sigma
              a -> a -> a
forall a. Subtractive a => a -> a -> a
- DiagonalMetric a -> [a] -> Int -> Int -> Int -> a
forall a.
Field a =>
DiagonalMetric a -> [a] -> Int -> Int -> Int -> a
gammaDiagonal DiagonalMetric a
m [a]
x Int
rho Int
nu Int
lam
              a -> a -> a
forall a. Multiplicative a => a -> a -> a
* DiagonalMetric a -> [a] -> Int -> Int -> Int -> a
forall a.
Field a =>
DiagonalMetric a -> [a] -> Int -> Int -> Int -> a
gammaDiagonal DiagonalMetric a
m [a]
x Int
lam Int
mu Int
sigma
          | Int
lam <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Subtractive a => a -> a -> a
- Int
1]
          ]
   in a
dMu a -> a -> a
forall a. Subtractive a => a -> a -> a
- a
dNu a -> a -> a
forall a. Additive a => a -> a -> a
+ a
quad

-- | Ricci @R_{σν} = R^ρ_{σρν}@.
ricciDiagonal ::
  (Field a, FromInteger a) =>
  DiagonalMetric a ->
  [a] ->
  Int ->
  Int ->
  a
ricciDiagonal :: forall a.
(Field a, FromInteger a) =>
DiagonalMetric a -> [a] -> Int -> Int -> a
ricciDiagonal DiagonalMetric a
m [a]
x Int
sigma Int
nu =
  let n :: Int
n = DiagonalMetric a -> Int
forall a. DiagonalMetric a -> Int
dmDim DiagonalMetric a
m
   in [a] -> a
forall a (f :: * -> *). (Additive a, Foldable f) => f a -> a
sum [DiagonalMetric a -> [a] -> Int -> Int -> Int -> Int -> a
forall a.
(Field a, FromInteger a) =>
DiagonalMetric a -> [a] -> Int -> Int -> Int -> Int -> a
riemannDiagonal DiagonalMetric a
m [a]
x Int
rho Int
sigma Int
rho Int
nu | Int
rho <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Subtractive a => a -> a -> a
- Int
1]]