A computational notebook on rotation geometry

Rotations2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "%config InlineBackend.figure_formats = ['svg']\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import init_session\n",
    "#init_session(quiet=True)\n",
    "from IPython.display import display, HTML\n",
    "\n",
    "# Hack: a function for centering equation output via HTML tags\n",
    "def center(eq1,eq2=None):\n",
    "    lxeq2 = \"\"\",\\qquad \"\"\"+latex(eq2) if eq2 is not None else ''\n",
    "    return HTML('<center>$'+latex(eq1)+lxeq2+'$</center>')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def prep3Dplot(size,zscale=1):\n",
    "    fig = plt.figure(figsize=(size,size))\n",
    "    axes = fig.add_subplot(111, projection='3d')\n",
    "    axes.set_position([0, 0, 1, 1])\n",
    "    axes.dist = 5\n",
    "    axes.set_axis_off()\n",
    "    axes.set_zlim3d(-zscale,zscale)\n",
    "    #axes.set_aspect('equal')\n",
    "    return fig, axes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mpl_toolkits.mplot3d import art3d\n",
    "from matplotlib.tri import Triangulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Spherical Projections"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "apply e.g. to component arrays of complex / quaternion\n",
    "Parallel or orthogonal projection simply drops the 1st component(s) of array.\n",
    "Central projection wrt 1st component(s) of array.\n",
    "Stereographic projection wrt 1st component(s) of array.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Spherical projections\n",
    "projections = {\n",
    "    \"parallel\":      lambda q: q[1:],\n",
    "    \"central\":       lambda q: q[1:]/q[0],\n",
    "    \"stereographic\": lambda q: q[1:]/(q[0]+1)\n",
    "}\n",
    "\n",
    "def SpherePrj(projection,q=None):\n",
    "    \"\"\"Spherical projection of ND vector q wrt 1st component.\"\"\"\n",
    "    P = projections[projection]\n",
    "    return P if q is None else P(q)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dotxy(N):\n",
    "    \"\"\"Return hemi-flat-torus-in-S3 as NxN meshgrid quaternions.\"\"\"\n",
    "    rangePi = np.arange(-N, N+1)*(np.pi/2/N)\n",
    "    ax,ay = np.meshgrid(rangePi,rangePi)\n",
    "    sax,cax = np.sin(ax), np.cos(ax)\n",
    "    say,cay = np.sin(ay), np.cos(ay)\n",
    "    # quaternion multiplication (cax,sax,0,0)*(cay,0,say,0)\n",
    "    return np.array([cax*cay, sax*cay, cax*say, sax*say])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fig,axes = prep3Dplot(8,1)\n",
    "axes.xaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))\n",
    "axes.yaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))\n",
    "axes.zaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))\n",
    "axes.dist = 6\n",
    "axes.elev = 45\n",
    "axes.azim = -65\n",
    "surf = axes.plot_surface(*SpherePrj(\"parallel\",dotxy(24)),\n",
    "                         rstride=2,cstride=2,\n",
    "                         cmap='viridis',alpha=0.8,\n",
    "                         edgecolor='black')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Rotation Representations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This  is a notebook narrative on the geometric and algebraic representation of rotations.\n",
    "\n",
    "Executable code is included for all results.\n",
    "\n",
    "The focus is on **$3D$** with an eye to **$ND$**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What's the point?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Points and their positions are the basis of Euclidean geometry. Newtonian dynamics identifies particles as point masses.\n",
    "\n",
    "In conceiving the field of **geometric mechanics**, Louis Poinsot considered rigid bodies as geometric objects in their own right rather than as collections of point particles."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">Everyone makes for himself a clear idea of the motion of a _point_, that is to say, of the motion of a corpuscle which one supposes to be infinitely small, and which one reduces by thought in some way to a mathematical point.\n",
    "\n",
    ">Tout le monde se fait une ideé claire du mouvement d'un _point_, c'est-à-dire du mouvement d'un corpuscle qu'on suppose infiniment petit, et qu'on réduit en quelque sorte par la pensée à un point mathématique.\n",
    "\n",
    ">—**Louis Poinsot**, Théorie nouvelle de la rotation des corps (1834)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    " By analogy, the whole can be considered a generalised point located in the space of possible positions. A geometric understanding of rotation "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Position, orientation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The position of a rigid body, its **'rigid position'**, consists of translational and rotational parts. This use of 'position' is more abstract than the position of a point in Euclidean space. The space of rigid positions is a non-Euclidean manifold, even if it is called the Special Euclidean group $SE(N)$. A rigid position can be treated as a point in this space. In 3D, rigid position is six dimensional. This idea of treating a rigid motion as the path of a point in a\n",
    "\n",
    "Rotational position is termed orientation or attitude. I tend to use orientation despite possible confusion with other meanings."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Non-linearity of position"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Position is non-linear; linearity requires addition but arbitrary addition of positions makes no sense. This is true even for point position, even if colloquially called 'linear' position, and even in 1D. \n",
    "\n",
    "Subtraction makes sense; the difference of two positions is the relative position from the subtrahend to the minuend (yes, I had to look that up). For point position, or just the translational position part of a rigid position, the difference **is** linear; the relative position vector. For orientation, the difference between a pair of attitudes is a rotation which **is not** linear."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Rotation as relative orientation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**R**elative **O**rien**TATION**, geddit?!??!?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3D Rotation 'arc-vectors'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although rotation is nonlinear, there is a tangible representation of 3D rotations that admits a visual vector-like addition operation.\n",
    "\n",
    "A 3D rotation can be visualised as an **arc-vector**; a directed arc of a great-circle on the sphere. It is vector-like in that it has direction and magnitude. It is like a free to move on its fixed great circle. The addition law for arc vectors modifies the parallelogram law of vector addition to a more general 'bow tie' arrangement."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Non-commutativity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The parallelogram law for adding two vectors starts with the vectors arranged tail-to-tail before they are translated, each along the other until, head-to-tail, they lie along the opposite two sides of the complete parallelogram. \n",
    "\n",
    "From that point, one vector is 'parallel-translated' along the other so that its tail coincides with the other's head. It doesn't matter which lies first and which second; the two possible orders create a parallelogram with the resultant tail-to-head vector-sum along the main diagonal the same either way.\n",
    "\n",
    "There is an alternative, more generalisable, way to view vector addition. The parallelogram law   putting focus on that starting poing. Focus instead on the head-to-tail meeting pointstart with \n",
    "\n",
    "Given a pair of arc-vectors, move each along its respective great circles such that they meet head-to-tail at one of the intersection points of the two great circles. Drawing a great circle arc from tail of the first to the head of the second gives the resultant arc-vector.\n",
    "\n",
    "The parallelogram law does not work directly for 'arc-vectors' on the sphere because spherical geometry is non-Euclidean. The commutative property is lost; the two orders of addition give different results. Then, although 'parallel-transport' can be interpreted on the sphere as rotating one arc-vector along the other, the resultant tail-to-head vector does not represent the compound rotation (it also\n",
    "\n",
    "\n",
    "\n",
    "in a vector-like fashion\n",
    "The translational part is \n",
    "\n",
    "It is not a vector in the linear algebra sense "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Tensor and rotor: matrix and quaternion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The two pillars of rotational mathematics are the **tensor** and **rotor** representations.\n",
    "hen resolved into components"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rotation in 3D is often visualised via a vector triad, $x,y,z$.\n",
    "They form a 'frame'; a mutually orthogonal set of basis vectors for 3D space.\n",
    "This generalises to ND, though visualisation is not easy for $N>=4$\n",
    "\n",
    "linear operator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import symbol,symbols,Symbol, Matrix,MatrixSymbol,MatMul,refine,Q,Eq,latex"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ND = Symbol('N')\n",
    "Rs = MatrixSymbol('\\\\mathbb{R}',ND,ND)\n",
    "RRTs = refine(Rs*Rs.T, Q.orthogonal(Rs))\n",
    "\n",
    "center(Eq(Rs*Rs.T, RRTs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Rs3 = Rs.subs(ND,3)\n",
    "RRTs3 = RRTs.subs(ND,3)\n",
    "\n",
    "Rs3e = Rs3.as_explicit()\n",
    "RRTs3e = RRTs3.as_explicit()\n",
    "\n",
    "Eq(MatMul(Rs3e,Rs3e.T), RRTs3e)\n",
    "#center(Eq(MatMul(Rs3e,Rs3e.T), RRTs3e))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "skew = lambda V: Matrix(3,3, lambda i,j: 0 if i==j else ((i-j+1)%3-1)*V[3-i-j])\n",
    "\n",
    "deskew = lambda M: Matrix(3,1, lambda i,j: M[(i-1)%3,(i-2)%3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Rs = MatrixSymbol('R',3,3)\n",
    "Rij = Matrix(3,3, lambda i,j: Symbol('r_{}'.format('xyz'[i]+'xyz'[j])))\n",
    "Rcontext = Q.orthogonal\n",
    "\n",
    "Sij = Matrix(3,3, lambda i,j: 0 if i==j else ((i-j+1)%3-1)*Symbol('s_{}'.format('XYZ'[3-i-j])))\n",
    "Aij = Matrix(3,3, lambda i,j: Symbol('c_{}'.format('xyz'[i])) if i==j else Symbol('a_{}'.format('xyz'[(3-i-j)])))\n",
    "\n",
    "Ss = (Rs-Rs.T)/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w,x,y,z = symbols('w,x,y,z')\n",
    "qs = MatrixSymbol('q',1,4)\n",
    "qi = Matrix(1,4,[w,x,y,z])\n",
    "center(Eq(Rs,Rij), Eq(qs,qi))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}