---
title: "Conformal geometry with Clifford algebra"
author: "Robin K. S. Hankin"
output: html_vignette
bibliography: clifford.bib
vignette: >
  %\VignetteIndexEntry{Conformal geometry with Clifford algebra}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---


```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library("clifford")
```

```{r out.width='20%', out.extra='style="float:right; padding:10px"',echo=FALSE}
knitr::include_graphics(system.file("help/figures/clifford.png", package = "clifford"))
```

\newcommand{\grad}[2]{\left\langle #1\right\rangle_{#2}}
\newcommand{\ak}{\grad{\mathbf{A}}{k}}
\newcommand{\ba}{\mathbf{a}}
\newcommand{\bb}{\mathbf{b}}
\newcommand{\bn}{\mathbf{n}}
\newcommand{\bA}{\mathbf{A}}
\newcommand{\bS}{\mathbf{S}}
\newcommand{\bX}{\mathbf{X}}
\newcommand{\bB}{\mathbf{B}}
\newcommand{\bx}{\mathbf{x}}
\newcommand{\by}{\mathbf{y}}

\newcommand{\einf}{\mathbf{e}_\infty}
\newcommand{\ezero}{\mathbf{e}_0}


To cite the `clifford` package in publications please use
@hankin2022_clifford.  This short document shows how conformal
geometry may be implemented using the `clifford` R package; it follows
@hildenbrand2013 and @perwass2009.  Here we work in
$\operatorname{Cl}(p,q)$, which Perwass denotes $\mathbb{G}_{p,q}$.
The set of grade-$k$ objects is denoted $\mathbb{G}_{p,q}^k$.  First
we define the IPNS and OPNS [inner product null space and outer
product null space] of $\mathbf{A}\in\mathbb{G}_{p,q}^k$.

$$
\begin{eqnarray}
\operatorname{IPNS}(\mathbf{A}) &=& \left\lbrace\bx\in\mathbb{G}^1_{p,q}\colon\bx\cdot\mathbf{A}  = 0\right\rbrace\\
\operatorname{OPNS}(\mathbf{A}) &=& \left\lbrace\bx\in\mathbb{G}^1_{p,q}\colon\bx\wedge\mathbf{A} = 0\right\rbrace
\end{eqnarray}
$$

Thus if $\mathbf{A}=\ak=\bigwedge_{i=1}^k\ba_i$, then
$\operatorname{OPNS}(\mathbf{A})=\operatorname{span}(\ba_1,\ldots,\ba_k)\subseteq\mathbb{G}_{p,q}$.
Further, $\operatorname{OPNS}(\mathbf{A})$ is linear in the sense that
$\bx,\by\in\operatorname{OPNS}(\mathbf{A})$ implies
$\alpha\bx+\beta\by\in\operatorname{OPNS}(\mathbf{A})$ for any real
numbers $\alpha,\beta$.  We will use this system to express a
geometrical object $\mathcal{G}\subset\mathbb{R}^3$ (such as a sphere
or a line) as a point $\mathbf{P}$ in $\operatorname{Cl}(4,1)$; the
idea is that the IPNS or OPNS of $\mathbf{P}$ is $\mathcal{G}$.  This
allows one to express geometrical ideas in pure Clifford formalism,
without needing to take a dangerous and unsightly basis.

To work in R, we set up some basic features of conformal geometry.
Specifically, we consider the three Euclidean basis vectors
$\mathbf{e}_1$, $\mathbf{e}_2$, $\mathbf{e}_3$ together with two
additional basis vectors $\mathbf{e}_+$ and $\mathbf{e}_-$ obeying
$\mathbf{e}_+^2=1$, $\mathbf{e}_-^2=0$,
$\mathbf{e}_+\cdot\mathbf{e}_-=0$.  If we define

\[
\ezero=\frac{1}{2}\left(\mathbf{e}_--\mathbf{e}_+\right),\qquad
\einf =\mathbf{e}_- +\mathbf{e}_+
\]

then we may say that $\ezero$ represents "the origin" and $\einf$
represents "the point at infinity".  We see that $\ezero^2=\einf^2=0$
and $\einf\cdot\ezero=-1$; the geometric product is given by
$\einf\ezero=-E-1$ and $\ezero\einf=E-1$ where $E=\einf\wedge\ezero$.
It is straightforward to implement these objects using the `clifford`
package:

```{r setdims}
dimension <- 3
options("maxdim" = dimension+2)  # paranoid safety measure
signature(dimension + 1,1)
eplus <- basis(dimension+1)
eminus <- basis(dimension + 2)

e0 <-  (eminus - eplus)/2
einf <- eminus + eplus
E <- e0 ^ einf
```

So

```{r showe0einfE}
e0
einf
E
```

# Points


With these definitions, we can consider Euclidean vectors
$\ba,\bb\in\mathbb{R}^3$ and conformal embeddings $\bA,\bB$ given by

\[
\bA=C(\ba)=\ba+\ba^2\einf/2 +\ezero\qquad\bB=C(\bb)=\bb+\bb^2\einf/2 +\ezero
\]

This is straightforward in package idiom:

```{r definepointfunc}
point <- function(x){ as.1vector(x) + sum(x^2)*einf/2 + e0 }
```

Thus `point()` takes an R vector of length 3 and returns its conformal
embedding.  For example, we may translate points $(1,2,5)$ and
$(2,2,2)$ to their conformal equivalent:

```{r pointexample}
a <- c(1,2,5)
b <- c(2,2,2)
point(a)
point(b)
```


It can be shown that $\bA\cdot\bB=-\left\lVert\ba-\bb\right\rVert^2/2$
[where the dot product is the Clifford inner product, `%.%`].  Package
idiom to verify this would be


```{r verifypointfunc}
c(conformal=drop(point(a) %.% point(b)), Euclidean = -sum((a-b)^2)/2)
```

showing that the two results match.

## Sphere, IPNS

We can define a sphere with center $\ba$ and radius $\rho$ as 

\[
\bS=C(\ba) - \rho^2\einf/2
\]

and then the sphere is just the inner product null space of $\bS$,
that is $\left\lbrace\bx\colon\bS\cdot\bx=0\right\rbrace$.  This is
straightforward to implement computationally.  Suppose we have a
sphere of radius 5, center $(1,2,3)$:

```{r definesphere}
sphere <- function(x,r){ point(x) - r^2*einf/2}
S <- sphere(1:3,5)  # center (1,2,3) radius 5:
S
```

Then object `S` is a conformal representation of such a sphere.  The
radius $\rho$ can be calculated from $\bS^2=\rho^2$.  Package idiom:

```{r verifyradius}
drop(S^2)   # 5^2 = 25
```
Finding the center of the sphere is slightly more involved;
Hildenbrand calls this the _sandwich_ product given by
$P=\bS\einf\bS$:

```{r calcsandichprod}
S*einf*S
```

Hildenbrand shows that the scaling factor is $-2$ so this gives us

```{r usescaling}
-S*einf*S/2
```

which we recognise as the point $(1,2,3)$:

```{r recog123}
point(1:3)
```


## Sphere, OPNS

We can also consider the OPNS for a sphere, defined in terms of four
points that lie on it.  For example, if our four points are $(0,0,0)$,
$(1,0,0)$, $(0,1,0)$, $(0,0,1)$ we would have a rather involved
algebraic calculation resulting in a sphere of radius
$\frac{\sqrt{3}}{2}$ and center
$\left(\frac{1}{2},\frac{1}{2},\frac{1}{2}\right)$.

The conformal representation for a sphere in OPNS would be 

\[
S^*=P_1\wedge P_2\wedge P_3\wedge P_4
\]

and points that lie on the sphere are the set

\[\left\lbrace\bx\colon\bx\wedge S^*=0\right\rbrace.\]

Observe again that this parameterization does not require one to take
a basis of $\mathbb{R}^3$, as all the results are presented in terms
of vector quantities: at no point does one need to consider components
or elements of any vector.  The R idiom for defining the sphere
touching $\left\lbrace P_1,P_2,P_3,P_4\right\rbrace$ is
straightforward:

```{r fourpoints}
origin <- point(c(0,0,0))
px <- point(c(1,0,0))
py <- point(c(0,1,0))
pz <- point(c(0,0,1))

(S <- origin ^ px ^ py ^ pz)
```

And this is the representation for a sphere (translating this into a
center and radius is not yet implemented).  Slightly slicker R idiom might be:

```{r definespherestar}
spherestar <- function(...){Reduce(`^`,list(...))}
spherestar(origin, px, py, pz)
```

As a verification, we may check that point
$\left(\frac{1}{2},\frac{1}{2},\frac{1}{2}+\frac{\sqrt{3}}{2}\right)$
is on the sphere as well:

```{r checksphere}
p <- point(c(1,1,1+sqrt(3))/2)
Mod(p ^ S)
```

Above we see that $p\wedge S$ is zero to within numerical error.
Again observe that this is established _without_ taking a basis of
$\mathbb{R}^3$.


# Planes

A plane is defined as $\hat{\bn} + d\einf$, where $\hat{\bn}$ is the
unit normal and $d$ the distance to the origin.

```{r defineplanefunc}
plane <- function(n,d){ as.1vector(n/sqrt(sum(n^2))) + d*einf}
```

For example, we consider the plane $\Pi$ with normal
$\mathbf{n}=\left(1,2,5\right)$ and distance 7.  We may use `plane()`
above but first need to calculate
$\hat{\mathbf{n}}=\mathbf{n}/\left|\mathbf{n}\right|$:

```{r generatend}
n <- c(1,2,5)
nhat <- n/sqrt(sum(n^2))
d <- 7
Pi <- plane(nhat,d)
Pi
```

Just to check, we may verify that some points known to lie in $\Pi$
are in fact in its IPNS.  We observe that the point
$7\hat{\mathbf{n}}\in\Pi$, and also that vectors $\mathbf{u}=(2,-1,0)$
and $\mathbf{v}=(5,0,-1)$ are orthogonal to the normal of $\Pi$.  So
$7\hat{\mathbf{n}} + \alpha\mathbf{u} + \beta\mathbf{v}\in\Pi$ for any
real numbers $\alpha,\beta$.  In R:

```{r generateuvp1p2p3}
u <- c(2,-1,0)
v <- c(5,0,-1)
P1 <- point(d*nhat)                
P2 <- point(d*nhat + 1.3*u + 3.44*v)
P3 <- point(d*nhat - 6.1*u + 1.02*v)
```

Above we use some made-up values of $\alpha,\beta$ to generate `P1`,
`P2`, `P3` which are known to lie in plane $\Pi$.  Then to verify that
these points lie in the IPNS of $\Pi$ we need to take the inner
product with the conformal representation of $\Pi$:

```{r verifyplaneipns}
c(drop(Pi %.% P1),drop(Pi %.% P2),drop(Pi %.% P3))
```

Above we see zero (to numerical precision), showing that we do indeed
have $7\hat{\bn} + \alpha\mathbf{u} +
\beta\mathbf{v}\in\operatorname{IPNS}(\Pi)$ for at least these values
of $\alpha$ and $\beta$.  Alternatively, a plane can be thought of as
a sphere that touches the point at infinity; thus the OPNS of plane is
given by

\[
\pi^\star = P_1\wedge P_2\wedge P_3\wedge\einf
\]

where $P_1,P_2,P_3$ are any points that lie in the plane.  To
illustrate this we will use the three points used in the IPNS above:

```{r plane3points}
Pi2 <- P1 ^ P2 ^ P3 ^ einf
```

Then for verification we can create another point known to be in the
plane and check that this is in the OPNS.  Below we will use `p4`
which is $d\hat{\mathbf{n}} + 7.6\mathbf{u} - 9.23\mathbf{v}$:

```{r checkonplane}
p4 <- point(d*nhat + 7.6*u - 9.23*v)
Mod(Pi2 ^ p4)
```

Above we see a very small result showing that `p4` is indeed in the
OPNS of plane `Pi2`.

# Circle

A circle is defined as the intersection of two spheres:

```{r definecirclefun}
circle <- function(S1,S2){  # IPNS
    S1 ^ S2
}
```

For example

```{r usecirclefun}
circle(sphere(1:3,5),sphere(c(1.1,2.1,3.4),6))
```

A circle may be represented in the OPNS by specifying three points
that lie on it:

\[ Z^*=P_1\wedge P_2\wedge P_3\]


```{r definecirclestar}
circlestar <- function(...){  # OPNS; A^B^C
	jj <- list(...)
    stopifnot(length(jj) == dimension)
    Reduce(`^`,lapply(jj,point))
}

(CIRC <- circlestar(c(1,2,3),c(5,6,3),c(8,8,-2)))
```

verify:


```{r usenumericstofindapointoncircle,echo=FALSE,cache=TRUE}
M <- rbind (c(1,2,3),c(5,6,3),c(8,8,-2))
badness_center <- function(pos){
  bad1 <-  sd(c(
  sum((pos-M[1,])^2),
  sum((pos-M[2,])^2),
  sum((pos-M[3,])^2)))   # pos should be equidistant from rows of M
  bad2 <- abs(det(sweep(M,2,pos))) # pos should be collinear with rows of M
  return(bad1+bad2)
  }

center <- nlm(badness_center,c(0,0,0))$estimate
badness_pointoncircle <- function(pos){ # pos is a point [that should be] on the circle
bad1 <-  var(c(sum((center-M[1,])^2), sum((center-M[2,])^2), sum((center-M[3,])^2),   sum((center-pos  )^2)))
bad2 <-  abs(det(sweep(M,2,pos))) # pos should be coplanar with M[i,]
return(bad1+bad2)
}

poc <- point(nlm(badness_pointoncircle,c(0,0,0))$estimate)
```

```{r verifycircleusingnum}
poc  # point on circle, found numerically [chunk omitted]
poc ^ CIRC
```

Above we see that `poc` is at least close to the circle from the small
magnitude of the terms in the wedge product.

# Lines and point pairs

A line is the intersection of two planes; in R:

```{r defline}
line <- function(P1,P2){ P1 ^ P2 }
```

and a "point pair" is the intersection of three spheres; in R:

```{r defpp}
pointpair <- function(S1,S2,S3){ S1 ^ S2 ^ S3 }
```

It is not at all obvious that three spheres intersect in a pair of
points; and still less obvious that the process is associative.
However, we may verify associativity explicitly:


```{r verifypointpairassoc}
S1 <- sphere(c(3,2,4),3)
S2 <- sphere(c(3,1,4),4)
S3 <- sphere(c(1,3,3),3)
(S1^S2)^S3
(S1^S2)^S3 == S1^(S2^S3)
```

```{r restoredefault, echo=FALSE}
options("maxdim" = NULL) # restore default
```


## References