Quaternions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### [Oregon Curriculum Network](http://www.4dsolutions.net/ocn) <br />\n",
    "[Discovering Math with Python](Introduction.ipynb)\n",
    "\n",
    "\n",
    "# Chapter 11: QUATERNIONS\n",
    "\n",
    "Quaternions were invented by Sir William Rowan Hamilton around 1843 and were considered a breakthrough.  \n",
    "\n",
    "In subsequent chapters Willard Gibbs and Oliver Heaviside came up with a vector concept that proved easier to use for many of the physics applications for which quaternions had originally been proposed.\n",
    "\n",
    "However quaternions and vectors together have become chief tools in accomplishing rotation, in computer graphics and games, robotics, rocketry.  We may have entered the realm of that proverbially most difficult of disciplines:  rocket science.\n",
    "\n",
    "Quaternions have some advantages over rotation matrices.  It's easier to slice up a rotation in a process called [SLERP](https://en.wikipedia.org/wiki/Slerp).\n",
    "\n",
    "You might think of quaternions as \"vectors on steroids\" if that helps.  They have some properties in common with vectors, and are conceived to have a \"vector part\" somewhat like complex numbers have a \"real part\".  However, like complex numbers, they're considered numbers in their own right, actually a superset of the complex, which in turn contain the reals and so on down to N, the counting numbers (natural numbers).\n",
    "\n",
    "Unit quaternions (w, x, y, z) with w\\*\\*2 + x\\*\\*2 + y\\*\\*2 + z\\*\\*2 == 1, form a group under multiplication.  Any two unit quaternions, when multiplied, produce a unit quaternion (Closure), and every unit quaternion has an inverse (w, -x, -y, -z), such that q \\* q\\*\\*-1 gives the unit quaternion (1, 0, 0, 0).\n",
    "\n",
    "Indeed the elements i, j, k from which quaternions are made, abet the i of complex number fame, with two more 2nd roots of -1.  All three, and their three inverses, engage in a kind of dance. \n",
    "\n",
    "Group elements are: {i, j, k, 1, -1, -i, -j, -k}.  Every product of two of these elements, is in this set (closure); associativity holds; every element has an inverse such that the two give a product of 1, and 1 serves as the neutral (identity) element.\n",
    "\n",
    "<a data-flickr-embed=\"true\"  href=\"https://www.flickr.com/photos/kirbyurner/36363327133/in/dateposted-public/\" title=\"Homework Assignment\"><img src=\"https://farm5.staticflickr.com/4353/36363327133_1962435626.jpg\" width=\"500\" height=\"375\" alt=\"Homework Assignment\"></a><script async src=\"//embedr.flickr.com/assets/client-code.js\" charset=\"utf-8\"></script>\n",
    "\n",
    "Just the identities we've been given: i\\*\\*2 = j\\*\\*2 = k\\*\\*2 = i \\* j \\* k = -1 are suffient to derive the above Cayley Table.  This table, in turn enables use to flesh out the \\_\\_mul\\_\\_ method, by taking sixteen products (w, x, y, z times each w', x', y', z') and substituting for products such as -i \\* k and k \\* j.  Some of the derivations for are table are shown above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What we'll plan to do, when rotating point vectors of polyhedron P, by n degrees around unit rotation vector q, is initialize the corresponding unit quaternion rQ using function rotation(n, q) and then multiply every point vector of P in a \"sandwich\" between rQ and its inverse:  \n",
    "\n",
    "new_Pv = rQ \\* Pv \\* ~rQ  (how we apply rotations to each position vector Pv).\n",
    "\n",
    "Remember our Polyhedron objects, as data structures, start with a set of faces pegged to vertexes, expressed as a dict of position Vectors.  We used both xyz and quadray notation to map the same set of points.\n",
    "\n",
    "Below is the source code for doing that.  Understanding multiplication involves developing the multiplication table for i, j, k which play the role of basis vectors in this language game. i\\*\\*2 == j\\*\\*2 == k\\*\\*2 == -1, and i \\* j \\* k == -1 as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "..\n",
      "----------------------------------------------------------------------\n",
      "Ran 2 tests in 0.008s\n",
      "\n",
      "OK\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<unittest.runner.TextTestResult run=2 errors=0 failures=0>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from math import cos, sin, radians, pi\n",
    "from qrays import Vector\n",
    "import unittest\n",
    "\n",
    "class Quaternion:\n",
    "    \n",
    "    def __init__(self, w, x, y, z):\n",
    "        \"\"\"\n",
    "        w is the scalar part;\n",
    "        x,y,z comprise a vector part as \n",
    "        the coefficients of i,j,k respectively\n",
    "        where i,j,k are the three basis vectors\n",
    "        described by Hamilton such that \n",
    "        i**2 == j**2 == k**2 == -1, and \n",
    "        i * j * k == -1 as well.\n",
    "        \"\"\"\n",
    "        self.w = w\n",
    "        self.x = x\n",
    "        self.y = y\n",
    "        self.z = z\n",
    "            \n",
    "    def __mul__(self, other):\n",
    "        \"\"\"\n",
    "        Derived by inter-multiplying all four terms in \n",
    "        each Quaternion to get 16 products and then \n",
    "        simplifying according to such rules as ij = k,\n",
    "        ji = -k.\n",
    "        \n",
    "        See:  https://youtu.be/jlskQDR8-bY Mathoma\n",
    "        'Quaternions explained briefly'\n",
    "        \"\"\"\n",
    "        a, b, c, d = self.w, self.x, self.y, self.z\n",
    "        e, f, g, h = other.w, other.x, other.y, other.z\n",
    "        w = (a*e - b*f - c*g - d*h)\n",
    "        x = (a*f + b*e + c*h - d*g)\n",
    "        y = (a*g - b*h + c*e + d*f)\n",
    "        z = (a*h + b*g - c*f + d*e)\n",
    "        return Quaternion(w, x, y, z)\n",
    "    \n",
    "    def __invert__(self):\n",
    "        return Quaternion( self.w,\n",
    "                          -self.x,\n",
    "                          -self.y,\n",
    "                          -self.z)\n",
    "        \n",
    "    def vector(self):\n",
    "        return Vector((self.x, self.y, self.z))\n",
    "        \n",
    "    def __eq__(self, other):\n",
    "        tolerance = 1e-8\n",
    "        return (abs(self.w - other.w) < tolerance and\n",
    "                abs(self.x - other.x) < tolerance and \n",
    "                abs(self.y - other.y) < tolerance and\n",
    "                abs(self.z - other.z) < tolerance)\n",
    "         \n",
    "    def __repr__(self):\n",
    "        return \"Quaternion({},{},{},{})\".format(self.w, \n",
    "                          self.x,\n",
    "                          self.y,\n",
    "                          self.z)\n",
    "\n",
    "def rotator(uV, α):\n",
    "    w = cos(α/2)\n",
    "    x = uV.x * sin(α/2)\n",
    "    y = uV.y * sin(α/2)\n",
    "    z = uV.z * sin(α/2)\n",
    "    return Quaternion(w, x, y, z)\n",
    "    \n",
    "class TestQuaternion(unittest.TestCase):\n",
    "    \n",
    "    def test_inverse(self):\n",
    "        Q = Quaternion(0, 0, -1, 0)\n",
    "        inverse = ~Q\n",
    "        self.assertEqual(Q * inverse, Quaternion(1,0,0,0))\n",
    "            \n",
    "    def test_x(self):\n",
    "        rV = Quaternion(0,1,0,0)\n",
    "        rQ = rotator(Vector((0,1,0)), pi)\n",
    "        newQ = rQ * rV * ~rQ\n",
    "        self.assertTrue(newQ.vector() == Vector((-1, 0, 0)))\n",
    "        \n",
    "    def test_360(self):\n",
    "        rV = Quaternion(0,1,0,0)\n",
    "        one_deg = radians(1)\n",
    "        rQ = rotator(Vector((0,1,0)), one_deg)\n",
    "        for _ in range(360):\n",
    "            rV = rQ * rV * ~rQ\n",
    "        self.assertTrue(rV.vector() == Vector((1, 0, 0)))\n",
    "        \n",
    "a = TestQuaternion()\n",
    "# unittest.TextTestRunner().run(a)\n",
    "\n",
    "suite = unittest.TestLoader().loadTestsFromModule(a)\n",
    "unittest.TextTestRunner().run(suite)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Back to Chapter 10: [Complex Numbers](Complex%20Numbers.ipynb) <br />\n",
    "Continue to Chapter 12: [The Mandelbrot Set](Mandelbrot%20Set.ipynb)<br />\n",
    "[Introduction / Table of Contents](Introduction.ipynb)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}