-- Hodge Laplacian in spherical coordinates

declare symbol r, θ, φ

def N : Integer := 3

def x : Vector MathExpr := [| r, θ, φ |]

def g_i_j : Matrix MathExpr := [| [| 1, 0, 0 |], [| 0, r^2, 0 |], [| 0, 0, r^2 * (sin θ)^2 |] |]_i_j
def g~i~j : Matrix MathExpr := [| [| 1, 0, 0 |], [| 0, 1 / r^2, 0 |], [| 0, 0, 1 / (r^2 * (sin θ)^2) |] |]~i~j

-- Exterior derivative
def d (A: Tensor MathExpr) : Tensor MathExpr := !(flip ∂/∂) x A

-- Hodge star operator
def hodge (A: DiffForm MathExpr) : DiffForm MathExpr :=
  let k := dfOrder A
   in withSymbols [i, j]
        sqrt (abs (M.det g_#_#)) *
        foldl
          (.)
          ((subrefs A (map 1#j_$1 (between 1 k))) . (subrefs (ε' N k) (map 1#i_$1 (between 1 N))))
          (map (\n -> g~(i_n)~(j_n)) (between 1 k))

-- Codifferential
def δ (A: DiffForm MathExpr) : DiffForm MathExpr :=
  let k := dfOrder A
   in ((-1)^(N * k + 1)) * hodge (d (hodge A))

-- Laplacian
def Δ (A: DiffForm MathExpr) : DiffForm MathExpr :=
  match dfOrder A as integer with
    | #0 -> δ (d A)
    | #N -> d (δ A)
    | _ -> d (δ A) + δ (d A)

-- Apply to scalar function
def f := function (r, θ, φ)

assertEqual "Laplacian in spherical coordinates"
  (Δ f)
  (∂/∂ (∂/∂ f r) r + 2 * (∂/∂ f r) / r + (∂/∂ (∂/∂ f θ) θ) / r^2 + cos θ * (∂/∂ f θ) / (r^2 * sin θ) + (∂/∂ (∂/∂ f φ) φ) / (r^2 * (sin θ)^2))