∑MathIntermediate

Riemannian Metrics & Geometry

Key Points

•
A Riemannian metric assigns an inner product to each tangent space, giving you a way to measure lengths and angles on curved spaces (manifolds).
•
Geodesics are the “straightest possible” paths under the metric; on a sphere, great circles are geodesics.
•
In coordinates, a metric is represented by a positive-definite matrix $g_{ij}$ (x), and lengths/angles use this matrix the way Euclidean geometry uses dot products.
•
Christoffel symbols, computed from $g_{ij}$ , describe how vectors change as you move and define the geodesic equations.
•
The gradient of a function on a manifold depends on the metric: grad $f = g^{- 1}$ ∇f, not just the usual coordinate gradient.
•
Volumes and integrals use the factor sqrt(det g); this adjusts for how coordinates stretch under the metric.
•
Numerical methods (finite differences, matrix inversion, and Runge–Kutta) let us compute geodesics, lengths, and gradients in C++.
•
Watch for coordinate singularities (like the poles on a sphere) and ensure your metric remains positive-definite for stable computations.

Prerequisites

→Linear algebra: inner products and positive-definite matrices — Metrics are inner products represented by symmetric positive-definite matrices g_{ij} in coordinates.
→Multivariable calculus: partial derivatives and chain rule — Christoffel symbols and gradients require derivatives of metric components and functions.
→Ordinary differential equations (ODEs) — Geodesics solve a second-order ODE that is integrated numerically.
→Coordinate systems and change of variables — Metrics transform under reparameterization; computations are chart-dependent.
→Numerical analysis basics — Finite differences, numerical integration, and stability are needed for practical computation.
→Vector calculus on surfaces (optional) — Understanding tangent vectors, surface parametrizations, and induced metrics helps with intuition.

Detailed Explanation

Tap terms for definitions

01Overview

Hook: Imagine hiking on Earth and wanting the shortest walking route between two cities. On a flat map, straight lines seem shortest, but on Earth (a sphere), the shortest route is a segment of a great circle. To capture such behavior mathematically, we need a way to measure distances and angles on curved spaces. Concept: A Riemannian metric equips each point of a manifold (a smoothly curved space) with a local ruler and protractor—an inner product on the tangent space. This inner product varies smoothly from point to point, letting us define lengths of curves, angles between directions, areas/volumes, gradients, curvature, and more. Example: On the unit sphere S^2 with spherical coordinates (θ, φ), the metric is ds^2 = dθ^2 + sin^2θ dφ^2. This tells us that a small change in longitude φ counts less near the poles (because circles of latitude are shorter), and more near the equator.

02Intuition & Analogies

Hook: Think of putting a tiny flat ruler on the surface at each point. Depending on where you are, the same coordinate step may correspond to different physical lengths. Concept: The Riemannian metric is a field of these tiny rulers—formally, a smoothly varying inner product on each tangent plane. It answers: How long is a tiny step? What’s the angle between two directions here? Analogy 1 (stretchy fabric): Picture a stretchy fabric with a grid printed on it. If the fabric is stretched more in one area and direction, a 1 cm grid square might become 2 cm by 1 cm; locally, distances in different directions scale differently. The metric is the information about this local stretching and shearing. Analogy 2 (GPS and maps): A map projection distorts areas and distances. The metric encodes how much a small patch on the map differs from the true patch on Earth, giving a local conversion rate. Example: On a sphere, going 1 degree east near the equator moves you a lot; near the pole, it’s barely any distance. The factor sinθ in ds^2 = dθ^2 + sin^2θ dφ^2 precisely captures that change.

03Formal Definition

A Riemannian metric g on a smooth n-dimensional manifold M assigns to each point p ∈ M a positive-definite, symmetric bilinear form

g_{p}

T_{p} M

T_{p} M

→ ℝ, varying smoothly with p. In local coordinates (

x^{1}

,…,

x^{n}

), g is represented by a symmetric positive-definite matrix field

g_{ij}

(x) so that for a tangent vector

v = v^{i}

∂/∂

x^{i}

, its squared length is g(v,v) =

g_{ij}

(x)

v^{i}

v^{j}

. The line element (infinitesimal squared distance) is

d s^{2} = g_{ij} (x

) d

x^{i}

x^{j}

. For a smooth curve

γ :

[

a, b] \to M

with velocity

γ^{'} (t)

\overset{x}{˙}^{i}

(t) ∂/∂

x^{i}

, its length is L(

γ)

= ∫_

a^{b}

√(

g_{ij}

(x(t))

\overset{x}{˙}^{i}

\overset{x}{˙}^{j}

) dt. The Levi–Civita connection ∇ is the unique torsion-free connection compatible with g (∇

g = 0

). In coordinates it is given by Christoffel symbols Γ^

k_{ij}

\frac{1}{2}

g^{k l}

(∂_i

g_{j l}

+ ∂_j

g_{i l}

− ∂_l

g_{ij}

). Geodesics satisfy

\overset{x}{¨}^{k}

+ Γ^

k_{ij}

\overset{x}{˙}^{i}

\overset{x}{˙}^{j}

= 0. The Riemann curvature tensor R encodes how much parallel transport depends on path and measures intrinsic curvature.

04When to Use

Measuring lengths/angles/areas on curved domains: Use a Riemannian metric when your space is not flat but you still need Euclidean-like notions locally (robot motion on curved surfaces, graphics on meshes, cartography).
Shortest paths (geodesics): Compute navigation routes on spheres or curved surfaces, intrinsic mesh processing, and shape analysis.
Optimization on manifolds: For problems constrained to lie on spheres, Stiefel/Grassmann manifolds, etc., the metric defines gradients, step sizes, and conjugacy for Riemannian optimization.
Physics and PDEs on manifolds: Heat flow, wave equations, and Laplace–Beltrami eigenproblems depend on g; simulations on surfaces use the metric for discretization.
Machine learning and statistics: Information geometry equips parameter spaces with Fisher information metrics to measure distances between distributions and do natural gradient methods.
Coordinate changes: When you reparameterize a surface, the pullback metric transforms correctly and preserves intrinsic quantities (length, area, curvature).

⚠️Common Mistakes

Treating coordinate differences as physical lengths: Always measure lengths with the metric ds^2 = g_{ij} dx^i dx^j, not naive Euclidean norms of coordinate changes.
Ignoring positive-definiteness: A valid Riemannian metric must be symmetric and positive-definite; numerical approximations that break this (e.g., due to round-off) cause unstable geodesic/gradient computations.
Forgetting coordinate singularities: Spherical coordinates break down at the poles (sinθ = 0). Avoid dividing by sinθ there; switch charts or regularize.
Mixing covariant/contravariant indices: Gradients in coordinates are covectors (∂_i f), while vectors use raised indices. Convert with g^{ij} to get grad f.
Not normalizing time when comparing curve lengths: Parameter speed affects the integrand, but length is parameterization-invariant only if you integrate √(g_{ij} \dot{x}^i \dot{x}^j) dt correctly.
Numerical differentiation too coarse/fine: Finite-difference steps that are too large lose accuracy; too small amplify floating-point errors. Use scale-aware step sizes and central differences.
Using explicit integrators with large steps: Geodesic ODEs can be stiff near singularities; prefer adaptive or sufficiently small steps (e.g., RK4) and monitor energy.

Key Formulas

Line Element

d s^{2} = g_{ij} (x) d x^{i} d x^{j}

Explanation: This is the infinitesimal squared distance in local coordinates. The matrix $g_{ij}$ determines how coordinate differentials turn into physical lengths.

Inner Product

⟨ v, w ⟩_{g} = g_{ij} (x) v^{i} w^{j}

Explanation: At a point, the metric acts like a dot product on tangent vectors. It measures angles and lengths of vectors in the tangent space.

Curve Length

L (γ) = \int_{a}^{b} g_{ij} (x (t)) \overset{x}{˙}^{i} (t) \overset{x}{˙}^{j} (t) d t

Explanation: The length of a curve is the integral of its speed, where speed is measured by the metric. This generalizes Euclidean arc length to curved spaces.

Energy Functional

E (γ) = \frac{1}{2} \int_{a}^{b} g_{ij} (x (t)) \overset{x}{˙}^{i} (t) \overset{x}{˙}^{j} (t) d t

Explanation: Geodesics minimize energy (and length among constant-speed curves). This variational view leads to the geodesic equation.

Christoffel Symbols

Γ_{ij}^{k} = \frac{1}{2} g^{k l} (\partial_{i} g_{j l} + \partial_{j} g_{i l} - \partial_{l} g_{ij})

Explanation: These coefficients encode how the coordinate basis changes and define the Levi–Civita connection. They are computed from g and its first derivatives.

Geodesic Equation

\overset{x}{¨}^{k} + Γ_{ij}^{k} \overset{x}{˙}^{i} \overset{x}{˙}^{j} = 0

Explanation: A geodesic has zero covariant acceleration. Solving this ODE gives locally shortest paths under the metric.

Riemannian Gradient

grad f = g^{ij} \partial_{j} f \partial_{i}

Explanation: To turn the coordinate gradient (a covector) into a vector, raise the index using the inverse metric $g^{ij}$ . This is the direction of steepest increase under g.

Laplace–Beltrami

Δ f = \frac{1}{∣ g ∣} \partial_{i} (∣ g ∣ g^{ij} \partial_{j} f)

Explanation: This generalizes the Laplacian to manifolds. It appears in diffusion, heat flow, and spectral geometry on curved spaces.

Volume Element

d V = ∣ g ∣ d x^{1} \dots d x^{n}

Explanation: Integrals over the manifold require the Jacobian factor √|g|, the square root of the determinant of the metric matrix.

Pullback Metric

g = J^{⊤} J

Explanation: For an embedding r(x): U $\subset$ $R^{n}$ $\to$ $R^{m}$ , the induced metric equals the Jacobian’s Gram matrix. It measures lengths as inherited from the ambient Euclidean space.

Sphere Metric

d s_{S^{2}}^{2} = d θ^{2} + sin^{2} θ d φ^{2}

Explanation: In spherical coordinates on the unit sphere, longitude steps shrink by sinθ near the poles, matching geometric intuition.

Complexity Analysis

Let n be the manifold’s dimension (

n = 2

for a surface chart) and let evaluating

g_{ij}

at a point cost O(

C_{g}

). For embedded surfaces r(u) in

R^{m}

with analytic partials, computing

g = J^{T}

J requires O(m

n^{2}

) arithmetic. If J is approximated via finite differences with k samples per coordinate, cost grows by O(k m n) per evaluation. Computing the inverse metric

g^{ij}

via a direct solver is O(

n^{3}

) in general, but for small fixed n (e.g.,

n = 2

) it is constant-time. Christoffel symbols require first derivatives ∂_k

g_{ij}

(O(

n^{2}

) quantities for each k). Using central finite differences, each derivative needs two metric evaluations, so evaluating all Γ^

k_{ij}

incurs O(

n^{4}

C_{g}

) work; with

n = 2

this is still constant but with a larger constant factor. Geodesic integration solves a 2n-dimensional first-order ODE using, say, Runge–Kutta 4 (RK4). Each step evaluates Γ and forms quadratic terms in velocities, costing O(

n^{3}

) dominated by symbol evaluation. With S integration steps, total time is O(S

n^{3}

+ S

n^{4}

C_{g}

) (effectively O(S) for small fixed n). Memory usage is O(

n^{2}

) to store g,

g^{- 1}

, and Γ at the current point, plus O(1) for the RK4 state; overall space is O(1) for fixed n. Curve length via numerical quadrature over M samples evaluates the integrand √(

g_{ij}

\overset{x}{˙}^{i}

\overset{x}{˙}^{j}

) at each sample: O(M

C_{g}

) time and O(1) space. Gradient computations use

g^{- 1}

times the coordinate gradient; that’s O(

n^{3}

) for inversion (constant for

n = 2

) and O(1) additional for multiplying by ∂f, so per-point cost is effectively constant in low dimensions. Stability-wise, errors accumulate linearly with step count for first-order methods and with higher-order rates for RK4; adaptive step sizing or re-normalization (e.g., constant-speed enforcement) helps maintain accuracy, especially near coordinate singularities where metric entries vary rapidly.

Code Examples

Compute metric, inner products, and curve length on the unit sphere S^2

1 #include <iostream>
2 #include <cmath>
3 #include <vector>
4 #include <array>
5 #include <iomanip>
6 
7 // This example uses spherical coordinates (theta, phi) on the unit sphere S^2.
8 // Metric: g = [[1, 0], [0, sin^2(theta)]]
9 // We compute:
10 //  - The metric tensor g at a point
11 //  - Inner product and norm of tangent vectors using g
12 //  - The length of a parametrized curve via numerical quadrature
13 
14 struct MetricSphere {
15     // Return metric matrix g at (theta, phi)
16     static std::array<std::array<double,2>,2> g(double theta, double /*phi*/) {
17         double s = std::sin(theta);
18         return {{{{1.0, 0.0}}, {{0.0, s*s}}}};
19     }
20 
21     // Inner product <v, w>_g at (theta, phi) for tangent vectors v, w in coordinate basis
22     static double inner(double theta, double phi, const std::array<double,2>& v, const std::array<double,2>& w) {
23         auto G = g(theta, phi);
24         // v^T G w = sum_{i,j} v^i g_{ij} w^j
25         return v[0]*(G[0][0]*w[0] + G[0][1]*w[1]) + v[1]*(G[1][0]*w[0] + G[1][1]*w[1]);
26     }
27 
28     // Norm induced by g
29     static double norm(double theta, double phi, const std::array<double,2>& v) {
30         return std::sqrt(inner(theta, phi, v, v));
31     }
32 };
33 
34 // Numerical integration (trapezoidal rule) of curve length
35 // gamma: [0,1] -> (theta, phi), provided as a function returning (theta(t), phi(t)) and derivative (theta'(t), phi'(t)).
36 struct Curve {
37     std::function<std::array<double,2>(double)> x;     // position (theta, phi)
38     std::function<std::array<double,2>(double)> dxdt;  // velocity (theta', phi')
39 };
40 
41 double curve_length(const Curve& gamma, int N) {
42     // Integrate L = ∫ sqrt( g_ij(x(t)) x'^i x'^j ) dt over t in [0,1]
43     double L = 0.0;
44     double dt = 1.0 / N;
45     for (int k = 0; k <= N; ++k) {
46         double t = k * dt;
47         auto pos = gamma.x(t);
48         auto vel = gamma.dxdt(t);
49         double speed = MetricSphere::norm(pos[0], pos[1], vel);
50         double w = (k == 0 || k == N) ? 0.5 : 1.0; // trapezoid weights
51         L += w * speed;
52     }
53     return L * dt;
54 }
55 
56 int main() {
57     std::cout.setf(std::ios::fixed); std::cout << std::setprecision(6);
58 
59     // Example 1: Metric and inner product at the equator (theta = pi/2)
60     double theta = M_PI / 2.0, phi = 0.0;
61     auto G = MetricSphere::g(theta, phi);
62     std::cout << "Metric g at (theta=pi/2, phi=0):\n";
63     std::cout << "[ [" << G[0][0] << ", " << G[0][1] << "],\n  [" << G[1][0] << ", " << G[1][1] << "] ]\n";
64 
65     // Tangent vectors in coordinates: v = (v_theta, v_phi)
66     std::array<double,2> v{1.0, 0.0}; // pointing in +theta direction
67     std::array<double,2> w{0.0, 1.0}; // pointing in +phi direction
68     std::cout << "<v,w>_g = " << MetricSphere::inner(theta, phi, v, w) << " (should be 0 at equator)\n";
69     std::cout << "||w||_g = " << MetricSphere::norm(theta, phi, w) << " (at equator equals 1)\n";
70 
71     // Example 2: Length of the equator curve gamma(t) = (theta=pi/2, phi=2*pi*t), expected length 2*pi
72     Curve equator;
73     equator.x = [](double t){ return std::array<double,2>{ M_PI/2.0, 2.0*M_PI*t }; };
74     equator.dxdt = [](double /*t*/){ return std::array<double,2>{ 0.0, 2.0*M_PI }; };
75 
76     double L = curve_length(equator, 2000);
77     std::cout << "Approx length of equator = " << L << " (expected ~ " << 2.0*M_PI << ")\n";
78 
79     return 0;
80 }
81

We use the analytic metric for S^2 to compute inner products and norms of tangent vectors. Then we integrate the speed along a curve using the metric to get its length. For the equator, the exact length is 2π, and the trapezoidal rule numerically approximates it closely with enough samples.

Time: Curve length with N samples is O(N); metric and inner product at each sample are O(1).Space: O(1) additional space beyond the output and a few scalars.

Integrate the geodesic equation on S^2 using RK4 and Christoffel symbols

1 #include <iostream>
2 #include <cmath>
3 #include <array>
4 #include <vector>
5 #include <iomanip>
6 
7 // Geodesic ODE on the unit sphere S^2 in spherical coordinates (theta, phi):
8 // x^1 = theta, x^2 = phi
9 // Nonzero Christoffel symbols:
10 //   Gamma^theta_{phi phi} = -sin(theta) * cos(theta)
11 //   Gamma^phi_{theta phi} = Gamma^phi_{phi theta} = cot(theta) = cos(theta)/sin(theta)
12 // ODE system in first-order form for state y = [theta, phi, vtheta, vphi]:
13 //   d/dt theta = vtheta
14 //   d/dt phi   = vphi
15 //   d/dt vtheta = -Gamma^theta_{ij} v^i v^j
16 //   d/dt vphi   = -Gamma^phi_{ij}   v^i v^j
17 
18 struct ChristoffelSphere {
19     static void gamma(double theta, double /*phi*/, double& Gt_pp, double& Gp_tp) {
20         double s = std::sin(theta), c = std::cos(theta);
21         Gt_pp = -s * c;           // Gamma^theta_{phi phi}
22         Gp_tp = (std::abs(s) < 1e-12) ? 0.0 : (c / s); // Gamma^phi_{theta phi} (handle near poles cautiously)
23     }
24 };
25 
26 struct State { double th, ph, vth, vph; };
27 
28 State deriv(const State& y) {
29     double Gt_pp, Gp_tp; ChristoffelSphere::gamma(y.th, y.ph, Gt_pp, Gp_tp);
30     // velocities: v^theta = vth, v^phi = vph
31     double dth = y.vth;
32     double dph = y.vph;
33     double dvth = - ( Gt_pp * y.vph * y.vph );
34     double dvph = - ( 2.0 * Gp_tp * y.vth * y.vph ); // since only Gamma^phi_{theta phi} and symmetric term are nonzero
35     return {dth, dph, dvth, dvph};
36 }
37 
38 State rk4_step(const State& y, double h) {
39     auto add = [](const State& a, const State& b, double s){ return State{a.th + s*b.th, a.ph + s*b.ph, a.vth + s*b.vth, a.vph + s*b.vph}; };
40     State k1 = deriv(y);
41     State k2 = deriv(add(y, k1, h*0.5));
42     State k3 = deriv(add(y, k2, h*0.5));
43     State k4 = deriv(add(y, k3, h));
44     State out;
45     out.th  = y.th  + (h/6.0)*(k1.th  + 2*k2.th  + 2*k3.th  + k4.th);
46     out.ph  = y.ph  + (h/6.0)*(k1.ph  + 2*k2.ph  + 2*k3.ph  + k4.ph);
47     out.vth = y.vth + (h/6.0)*(k1.vth + 2*k2.vth + 2*k3.vth + k4.vth);
48     out.vph = y.vph + (h/6.0)*(k1.vph + 2*k2.vph + 2*k3.vph + k4.vph);
49     // keep phi within [-pi, pi] for readability
50     if (out.ph > M_PI) out.ph -= 2.0*M_PI; else if (out.ph < -M_PI) out.ph += 2.0*M_PI;
51     return out;
52 }
53 
54 int main(){
55     std::cout.setf(std::ios::fixed); std::cout << std::setprecision(6);
56 
57     // Start at equator (theta = pi/2, phi = 0) with velocity tangent along +phi.
58     // The equator is a geodesic, so theta should remain ~ pi/2.
59     State y{M_PI/2.0, 0.0, 0.0, 1.0}; // unit speed in phi at equator
60 
61     double h = 1e-3; // time step
62     int steps = 4000;
63 
64     for (int i = 0; i <= steps; ++i) {
65         if (i % 500 == 0) {
66             std::cout << "t=" << (i*h)
67                       << "  theta=" << y.th
68                       << "  phi=" << y.ph
69                       << "  vtheta=" << y.vth
70                       << "  vphi=" << y.vph << "\n";
71         }
72         y = rk4_step(y, h);
73     }
74 
75     return 0;
76 }
77

We hard-code the Christoffel symbols for S^2 and integrate the geodesic ODE using RK4. Starting along the equator with purely longitudinal velocity, the solution remains near θ = π/2, verifying the equator is a geodesic. The method generalizes by replacing the hard-coded Γ with values computed from an arbitrary metric.

Time: O(S) for S RK4 steps; each step is O(1) for S^2 with closed-form Γ.Space: O(1) to store the current state and temporary RK4 vectors.

Compute the Riemannian gradient on S^2 in spherical coordinates

1 #include <iostream>
2 #include <cmath>
3 #include <iomanip>
4 
5 // Riemannian gradient on S^2 (theta, phi) for a scalar function f(theta, phi):
6 // g^{-1} = diag(1, 1/sin^2(theta))
7 // grad f = (df/dtheta) * d/dtheta + (1/sin^2(theta)) * (df/dphi) * d/dphi
8 
9 struct GradSphere {
10     static std::pair<double,double> grad(double theta, double phi, 
11                                          double dfdtheta, double dfdphi) {
12         double s = std::sin(theta);
13         double inv_s2 = (std::abs(s) < 1e-12) ? 0.0 : 1.0/(s*s); // caution near poles
14         double gth = dfdtheta;
15         double gph = inv_s2 * dfdphi;
16         return {gth, gph}; // components in (theta, phi) basis
17     }
18 };
19 
20 // Example function: f(theta, phi) = cos(theta) (height relative to z-axis)
21 // Maximum at theta=0 (north pole), minimum at theta=pi (south pole)
22 
23 int main(){
24     std::cout.setf(std::ios::fixed); std::cout << std::setprecision(6);
25 
26     auto f = [](double th, double /*ph*/){ return std::cos(th); };
27     auto df = [](double th, double /*ph*/){
28         double dth = -std::sin(th); // df/dtheta
29         double dph = 0.0;           // df/dphi
30         return std::pair<double,double>(dth, dph);
31     };
32 
33     double th = 1.2;   // initial theta (~68.8 degrees)
34     double ph = 0.5;   // initial phi
35     double eta = 0.1;  // step size for gradient ascent
36 
37     for (int k = 0; k < 20; ++k) {
38         auto [dth, dph] = df(th, ph);
39         auto [gth, gph] = GradSphere::grad(th, ph, dth, dph);
40         // Gradient ascent step in coordinates
41         th += eta * gth;
42         ph += eta * gph;
43         // wrap phi
44         if (ph > M_PI) ph -= 2.0*M_PI; else if (ph < -M_PI) ph += 2.0*M_PI;
45         // clamp theta to (0, pi)
46         if (th < 1e-6) th = 1e-6; if (th > M_PI - 1e-6) th = M_PI - 1e-6;
47         std::cout << "iter=" << k << "  theta=" << th << "  phi=" << ph
48                   << "  f=" << f(th, ph) << "\n";
49     }
50 
51     return 0;
52 }
53

We compute the Riemannian gradient using the inverse metric. For f(θ,φ) = cos θ, the gradient points toward decreasing θ, i.e., toward the north pole when doing ascent on −f or descent on f. The example performs simple gradient ascent steps to increase f toward its maximum at θ = 0.

Time: O(T) for T iterations; each iteration is O(1) for evaluating derivatives and applying g^{-1}.Space: O(1) additional space.

1	#include <iostream>
2	#include <cmath>
3	#include <vector>
4	#include <array>
5	#include <iomanip>
6
7	// This example uses spherical coordinates (theta, phi) on the unit sphere S^2.
8	// Metric: g = [[1, 0], [0, sin^2(theta)]]
9	// We compute:
10	// - The metric tensor g at a point
11	// - Inner product and norm of tangent vectors using g
12	// - The length of a parametrized curve via numerical quadrature
13
14	struct MetricSphere {
15	// Return metric matrix g at (theta, phi)
16	static std::array<std::array<double,2>,2> g(double theta, double /phi/) {
17	double s = std::sin(theta);
18	return {{{{1.0, 0.0}}, {{0.0, s*s}}}};
19	}
20
21	// Inner product <v, w>_g at (theta, phi) for tangent vectors v, w in coordinate basis
22	static double inner(double theta, double phi, const std::array<double,2>& v, const std::array<double,2>& w) {
23	auto G = g(theta, phi);
24	// v^T G w = sum_{i,j} v^i g_{ij} w^j
25	return v[0](G[0][0]w[0] + G[0][1]w[1]) + v[1](G[1][0]w[0] + G[1][1]w[1]);
26	}
27
28	// Norm induced by g
29	static double norm(double theta, double phi, const std::array<double,2>& v) {
30	return std::sqrt(inner(theta, phi, v, v));
31	}
32	};
33
34	// Numerical integration (trapezoidal rule) of curve length
35	// gamma: [0,1] -> (theta, phi), provided as a function returning (theta(t), phi(t)) and derivative (theta'(t), phi'(t)).
36	struct Curve {
37	std::function<std::array<double,2>(double)> x; // position (theta, phi)
38	std::function<std::array<double,2>(double)> dxdt; // velocity (theta', phi')
39	};
40
41	double curve_length(const Curve& gamma, int N) {
42	// Integrate L = ∫ sqrt( g_ij(x(t)) x'^i x'^j ) dt over t in [0,1]
43	double L = 0.0;
44	double dt = 1.0 / N;
45	for (int k = 0; k <= N; ++k) {
46	double t = k * dt;
47	auto pos = gamma.x(t);
48	auto vel = gamma.dxdt(t);
49	double speed = MetricSphere::norm(pos[0], pos[1], vel);
50	double w = (k == 0 \|\| k == N) ? 0.5 : 1.0; // trapezoid weights
51	L += w * speed;
52	}
53	return L * dt;
54	}
55
56	int main() {
57	std::cout.setf(std::ios::fixed); std::cout << std::setprecision(6);
58
59	// Example 1: Metric and inner product at the equator (theta = pi/2)
60	double theta = M_PI / 2.0, phi = 0.0;
61	auto G = MetricSphere::g(theta, phi);
62	std::cout << "Metric g at (theta=pi/2, phi=0):\n";
63	std::cout << "[ [" << G[0][0] << ", " << G[0][1] << "],\n [" << G[1][0] << ", " << G[1][1] << "] ]\n";
64
65	// Tangent vectors in coordinates: v = (v_theta, v_phi)
66	std::array<double,2> v{1.0, 0.0}; // pointing in +theta direction
67	std::array<double,2> w{0.0, 1.0}; // pointing in +phi direction
68	std::cout << "<v,w>_g = " << MetricSphere::inner(theta, phi, v, w) << " (should be 0 at equator)\n";
69	std::cout << "\|\|w\|\|_g = " << MetricSphere::norm(theta, phi, w) << " (at equator equals 1)\n";
70
71	// Example 2: Length of the equator curve gamma(t) = (theta=pi/2, phi=2pit), expected length 2*pi
72	Curve equator;
73	equator.x = [](double t){ return std::array<double,2>{ M_PI/2.0, 2.0M_PIt }; };
74	equator.dxdt = [](double /t/){ return std::array<double,2>{ 0.0, 2.0*M_PI }; };
75
76	double L = curve_length(equator, 2000);
77	std::cout << "Approx length of equator = " << L << " (expected ~ " << 2.0*M_PI << ")\n";
78
79	return 0;
80	}
81