# Differential geometry

## Directional derivative

```use math.la: vector

function diffvh(h)
return |v,f,x| (f(x+h*v)-f(x-h*v))/(2*h)
end

diffv = diffvh(0.001)

# The standard basis
e0 = vector(1,0)
e1 = vector(0,1)

# A scalar field
f = |[x,y]| x*y

# Value and partial derivatives at (x=2,y=4)
p = vector(2,4)
print(f(p))
print(diffv(e0,f,p))
print(diffv(e1,f,p))

# A vector field
F = |[x,y]| vector(x*y,x+y)

# Value and partial derivatives at (x=2,y=4)
p = vector(2,4)
print(F(p))
print(diffv(e0,F,p))
print(diffv(e1,F,p))

# Tangential plane of f at p
T = |f,p| |x| f(p)+diffv(x-p,f,p)

p = vector(2,4)
g = T(f,p)
G = T(F,p)

# Note that diffv may be used as Gâteaux differential.
# One just needs to define the arithmetic operations
# for functions:
Function.sub = |f;g| |x| f(x)-g(x)
Function.rmul= |r;f| |x| r*f(x)
Function.div = |f;r| |x| f(x)/r
```

```use math.la: unit, vector

# Nabla operator
function nabla(f,p)
n = p.shape
E = list(0..n-1).map(|k| unit(n,k))
return vector(*E.map(|e| diffv(e,f,p)))
end

# Tangential plane of f at p
function T(f,p)
a = nabla(f,p)
y = f(p)
return |x| y+a*(x-p)
end

f = |[x,y]| x*y
g = T(f,vector(2,4))

# Taking the gradient of a vector valued function
# is possible, because we use the polymorphic array
# data type and thus can have vectors of vectors
# (i.e. jagged arrays).
F = |[x,y]| vector(x*y,x+y)
G = T(F,vector(2,4))
```

## Jacobian matrix

```use math.la: unit, vector, matrix

# The Jacobian matrix is the transposed gradient
function Jacobian(F,p)
return matrix(*nabla(F,p).list().map(|v| v.list())).T
end

# After simplification we obtain
function Jacobian(F,p)
n = p.shape
E = list(0..n-1).map(|k| unit(n,k))
return matrix(*E.map(|e| diffv(e,F,p).list())).T
end

# Tangential plane of F at p
function T(F,p)
J = Jacobian(F,p)
y = F(p)
return |x| y+J*(x-p)
end
```

## Volume

```use math.na: integral
use math.la.inversion: det
use math: pi, sin, cos

function volume(F,[x0,x1],[y0,y1],[z0,z1])
return integral(x0,x1,|x| integral(y0,y1,|y| integral(z0,z1,|z|
abs(det(Jacobian(F,vector(x,y,z))))
)))
end

torus = |{R}| |[r,p,t]| vector(
(R+r*cos(p))*cos(t),
(R+r*cos(p))*sin(t),
r*sin(p)
)

F = torus(R=2)
V = volume(F,[0,1],[0,2*pi],[0,2*pi])
print(V)

# R=2; r=1
# V = 2*pi^2*r^2*R
```

## Area

```use math.na: integral
use math.la.inversion: det
use math: sqrt, pi, sin, cos

function metric_tensor(F)
return fn|p|
J = Jacobian(F,p)
return J.T*J
end
end

function area(F,[x0,x1],[y0,y1])
g = metric_tensor(F)
return integral(x0,x1,|x| integral(y0,y1,|y|
sqrt(det(g(vector(x,y))))
))
end

torus = |{R,r}| |[p,t]| vector(
(R+r*cos(p))*cos(t),
(R+r*cos(p))*sin(t),
r*sin(p)
)

F = torus(R=2,r=1)
A = area(F,[0,2*pi],[0,2*pi])
print(A)

# R=2; r=1
# A = 4*pi^2*r*R
```