wball/Circle-Packings
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions
Circle-Packings/.ipynb_checkpoints/Apollonian Circle Packings-checkpoint.ipynb

689 lines
80 KiB
Text
Raw Normal View History

tried some new stuff
2021-03-08 07:43:47 -08:00
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Calculations to find the quadratic form of the octahedron packing"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$P$ is the matrix of the quadratic form corresponding to $h_1^2 + h_2^2 - \\tilde{b}b = 1$."
]
},
{
"cell_type": "code",
"execution_count": 49,
"metadata": {},
"outputs": [],
"source": [
"P = matrix([\n",
" [0, -1/2, 0, 0],\n",
" [-1/2, 0, 0, 0],\n",
" [0, 0, 1, 0],\n",
" [0, 0, 0, 1],\n",
"])"
]
},
{
"cell_type": "code",
"execution_count": 50,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[ 0 -2 0 0]\n",
"[-2 0 0 0]\n",
"[ 0 0 1 0]\n",
"[ 0 0 0 1]"
]
},
"execution_count": 50,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"P.inverse()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For an Apollonian packing, you can make a matrix, $W$, where the rows are the coordinates of each of the circles in a quadruple. In an octahedral packing, this is impossible, as there are six circles in a unit. You can try creating a $6\\times 4$ matrix, but, later, we end up needing to invert $WPW^T$, the product of a $6\\times 4$, $4\\times 4$, and $4\\times 6$ matrix, which is singular. Fortunately, the sextuples come in three pairs of circles. Each pair consists of circles that aren't tangent to each other. For whatever reason (we aren't sure yet why, but the Guettler and Mallows says this is the case), the average of the coordinates of the circles in a pair is the same across the sextuple. So, we can make a matrix where the first three rows are the coordinates of circles from different pairs, i.e. three mutually tangent circles, and the fourth is the average of the coordinates in a pair. From this, you can recover the coordinates for all the circles in a sextuple.\n",
"\n",
"Here $W$ is such a matrix computed for the $(0, 0, 1, 1, 2, 2)$ root sextuple."
]
},
{
"cell_type": "code",
"execution_count": 51,
"metadata": {},
"outputs": [],
"source": [
"W = matrix([\n",
" [2, 0, 0, 1],\n",
" [2, 0, 0, -1],\n",
" [-1, 1, 0, 0],\n",
" [3, 1, sqrt(2), 0]\n",
"])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here we have $$\n",
" WPW^T = \\left(\\begin{matrix}\n",
" 1 & -1 & -1 & -1\\\\\n",
" -1 & 1 & -1 & -1\\\\\n",
" -1 & -1 & 1 & -1\\\\\n",
" -1 & -1 & -1 & -1\n",
" \\end{matrix}\\right) = M\n",
"$$\n",
"So, inverting both sides, we get $$\n",
" (WPW^T)^{-1} = M^{-1}\n",
"$$ $$\n",
" (W^T)^{-1}P^{-1}W^{-1} = M^{-1}\n",
"$$ $$\n",
" P^{-1} = W^TM^{-1}W\n",
".$$\n",
"\n",
"Like in the case with Apollonian packings, this is true for any sextuple $W$. So, we can substitute an arbitrary sextuple for $W$ and it must be equal to $P^{-1}$, letting us derive some useful quadratic forms."
]
},
{
"cell_type": "code",
"execution_count": 52,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/plain": [
"[ 1 -1 -1 -1]\n",
"[-1 1 -1 -1]\n",
"[-1 -1 1 -1]\n",
"[-1 -1 -1 -1]"
]
},
"execution_count": 52,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"M = W * P * W.transpose()\n",
"M"
]
},
{
"cell_type": "code",
"execution_count": 53,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/plain": [
"[ 1/2 0 0 -1/2]\n",
"[ 0 1/2 0 -1/2]\n",
"[ 0 0 1/2 -1/2]\n",
"[-1/2 -1/2 -1/2 1/2]"
]
},
"execution_count": 53,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"M.inverse()"
]
},
{
"cell_type": "code",
"execution_count": 54,
"metadata": {},
"outputs": [],
"source": [
"bt1 = var('bt1')\n",
"b1 = var('b1')\n",
"h11 = var('h11')\n",
"h12 = var('h12')\n",
"bt2 = var('bt2')\n",
"b2 = var('b2')\n",
"h21 = var('h21')\n",
"h22 = var('h22')\n",
"bt3 = var('bt3')\n",
"b3 = var('b3')\n",
"h31 = var('h31')\n",
"h32 = var('h32')\n",
"b5_avg = var('b5_avg')\n",
"b_avg = var('b_avg')\n",
"h1_avg = var('h1_avg')\n",
"h2_avg = var('h2_avg')\n",
"\n",
"\n",
"W2 = matrix([\n",
" [bt1, b1, h11, h12],\n",
" [bt2, b2, h21, h22],\n",
" [bt3, b3, h31, h32],\n",
" [b5_avg, b_avg, h1_avg, h2_avg],\n",
"])"
]
},
{
"cell_type": "code",
"execution_count": 55,
"metadata": {},
"outputs": [],
"source": [
"D = W2.transpose() * M.inverse() * W2"
]
},
{
"cell_type": "code",
"execution_count": 56,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[ 0 -2 0 0]\n",
"[-2 0 0 0]\n",
"[ 0 0 1 0]\n",
"[ 0 0 0 1]"
]
},
"execution_count": 56,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"P.inverse()"
]
},
{
"cell_type": "code",
"execution_count": 57,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[b_avg == b1 + b2 + b3 - sqrt(2*b1*b2 + 2*(b1 + b2)*b3), b_avg == b1 + b2 + b3 + sqrt(2*b1*b2 + 2*(b1 + b2)*b3)]"
]
},
"execution_count": 57,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"quad = 2 * factor(simplify(D[1][1]))\n",
"solve(quad == 0, b_avg)"
]
},
{
"cell_type": "code",
"execution_count": 58,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"b1^2 + b2^2 + b3^2 - 2*b1*b_avg - 2*b2*b_avg - 2*b3*b_avg + b_avg^2"
]
},
"execution_count": 58,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"quad"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"So, we end up deriving the quadratic form $$\n",
" b_1^2 + b_2^2 + b_3^2 + b_{\\text{avg}}^2 - 2b_{\\text{avg}}(b_1 + b_2 + b_3) = 0\n",
".$$\n",
"\n",
"This means that, given three mutually tangent circles with curvatures $b_1,b_2,b_3$, there are two solutions for $b_{\\text{avg}}$, allowing us to derive two new sets of three mutually tangent circles with curvatures $b_1' = 2b_{\\text{avg}} - b_1$ etc."
]
},
{
"cell_type": "code",
"execution_count": 59,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1/2*b1*h11 - 1/2*b_avg*h11 - 1/2*b1*h1_avg - 1/2*b2*h1_avg - 1/2*b3*h1_avg + 1/2*b_avg*h1_avg + 1/2*b2*h21 - 1/2*b_avg*h21 + 1/2*b3*h31 - 1/2*b_avg*h31"
]
},
"execution_count": 59,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"factor(simplify(D[2][1]))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Calculations to find the generators for the \"Apollonian\" Group for the octahedral packing\n",
"\n",
"We can try the reflection group, since that is one way to derive the generators for the Apollonian group and see if it works. My worry is that it won't since the coordinate system is so different."
]
},
{
"cell_type": "code",
"execution_count": 60,
"metadata": {},
"outputs": [],
"source": [
"S1 = matrix([\n",
" [-1, 2, 2, 2],\n",
" [0, 1, 0, 0],\n",
" [0, 0, 1, 0],\n",
" [0, 0, 0, 1],\n",
"])\n",
"\n",
"S2 = matrix([\n",
" [1, 0, 0, 0],\n",
" [2, -1, 2, 2],\n",
" [0, 0, 1, 0],\n",
" [0, 0, 0, 1]\n",
"])\n",
"\n",
"S3 = matrix([\n",
" [1, 0, 0, 0],\n",
" [0, 1, 0, 0],\n",
" [2, 2, -1, 2],\n",
" [0, 0, 0, 1]\n",
"])\n",
"\n",
"# this is the only one that differs from the generators of the Apollonian group\n",
"S4 = matrix([\n",
" [1, 0, 0, 0],\n",
" [0, 1, 0, 0],\n",
" [0, 0, 1, 0],\n",
" [2, 2, 2, 3]\n",
"])"
]
},
{
"cell_type": "code",
"execution_count": 76,
"metadata": {},
"outputs": [],
"source": [
"W = matrix([\n",
" [2, 0, 0, 1],\n",
" [2, 0, 0, -1],\n",
" [-1, 1, 0, 0],\n",
" [3, 1, sqrt(2), 0]\n",
"])"
]
},
{
"cell_type": "code",
"execution_count": 85,
"metadata": {},
"outputs": [],
"source": [
"def draw_circles_from_coords(coords, color=(0,0,1)):\n",
" drawing = None\n",
" averages = coords[-1]\n",
" avgb = averages[1]\n",
" for row in coords[:-1]:\n",
" b = row[1]\n",
" if not b == 0:\n",
" x = row[2] / b\n",
" y = row[3] / b\n",
" if drawing is not None:\n",
" drawing += circle((x, y), 1 / b, rgbcolor=color)\n",
" else:\n",
" drawing = circle((x, y), 1 / b, rgbcolor=color)\n",
" newb = 2 * avgb - b\n",
" if not newb == 0:\n",
" newx = (2 * averages[2] - row[2]) / newb\n",
" newy = (2 * averages[3] - row[3]) / newb\n",
" if drawing is not None:\n",
" drawing += circle((newx, newy), 1 / newb, rgbcolor=color)\n",
" else:\n",
" drawing = circle((newx, newy), 1 / newb, rgbcolor=color)\n",
" return drawing"
]
},
{
"cell_type": "code",
"execution_count": 95,
"metadata": {},
"outputs": [],
"source": [
"def draw_circles_from_coords(coords, color=(0,0,1)):\n",
" drawing = None\n",
" for row in coords:\n",
" b = row[1]\n",
" if not b == 0:\n",
" x = row[2] / b\n",
" y = row[3] / b\n",
" if drawing is not None:\n",
" drawing += circle((x, y), 1 / b, rgbcolor=color)\n",
" else:\n",
" drawing = circle((x, y), 1 / b, rgbcolor=color)\n",
" return drawing"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"root = matrix([\n",
" [2, 0, 0, 1],\n",
" [2, 0, 0, -1],\n",
" [-1, 1, 0, 0],\n",
" [3, 1, 2, 0]\n",
"])\n",
"\n",
"S4 = matrix([\n",
" [1, 0, 0, 0],\n",
" [0, 1, 0, 0],\n",
" [0, 0, 1, 0],\n",
" [2, 2, 2, -1],\n",
"])\n",
"\n",
"drawing = draw_circles_from_coords(root)\n",
"branches = [root]\n",
"for i in range(10):\n",
" new_branches = []\n",
" for tree in branches:\n",
" drawing += draw_circles_from_coords(tree)\n",
" new_branches.append(S1 * tree)\n",
" new_branches.append(S2 * tree)\n",
" new_branches.append(S3 * tree)\n",
" new_branches.append(S4 * tree)\n",
" branches = new_branches\n",
"drawing"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Hmm. Doesn't look like it worked. And my crude code to draw the circles doesn't seem to be the problem either."
]
},
{
"cell_type": "code",
"execution_count": 93,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAkwAAAEECAYAAADah2tfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAA9hAAAPYQGoP6dpAABbr0lEQVR4nO3dd3iT1RcH8G+6aaFl0zJlbxkFBVkiywEFBWXJENmgoqgMRRBRVByIFGQoQ0D4KVuUPUQ2BWSjbBDKbksLXcn9/fE1tIWWpm3evHmT83mePNE2eXNo1nnPvfdck1IKQgghhBAiYx56ByCEEEII4ewkYRJCCCGEyIQkTEIIIYQQmZCESQghhBAiE5IwCSGEEEJkQhImIYQQQohMSMIkhBBCCJEJSZiEEEIIITIhCZN4gIkCTSaTSe9YhBBCCGfglYXbSktwNxEdHY2goCBER0frHYoQQgjhCJkWCKTC5OT++OMPtGnTBkWLFoXJZMKyZcsyvc/mzZtRu3Zt+Pr6oly5cpg9e7bmcQohhBCuTBImJxcXF4caNWogPDzcptufOXMGzz33HJo2bYoDBw5gyJAh6N27N9asWaNxpEIIIYTrMmVh810ZktOZyWTC0qVL0a5duwxvM2zYMKxatQqHDx++97NOnTohKioKq1evtulxYmJi7g3JBQYG5jRsIYQQwtllOiSXlTlMwgB27NiB5s2bp/lZq1at8Prr72LvXuDoUeDECeDsWeDiReDaNSAx0YLkZAWz2QSzGUhODgBwGRUr5oa3N+DhAXh6Al5eQFAQULQoUKoUUK4cUKkSUKMGULiwLv9cYUB37/L1d/lyyuXSJeDqVSAhAUhKAsxmvt68vYGAACAkhJeiRVOuS5Tg61KIzCgF3LyZ8lpLfR0bCyQn85KUxGuzOeUzz/o69PICcudO/7WYPz8gS2Rcn1SYDORhFSaLBdi+HWjTZiby5m0JH5+SuHWLHwbx8Qp8mtO+oz08AF9fwGxOQGLiXQBmAJb/LgoFChSGUh6wWHh8sxlITOSHSnp8fAB/fyZVBQsymXr6aSAsDJBClXu6cwf46y8gIiLlcvQoX0tWgYH84ilSBPDzw70k3foldvt2SmKVmJhyv4AAoFYtIDQUqFOH1xUqSBLlrm7dAvbtA/bvT0nIrUnR/a8dgElOSAhff/cnRp6e+O/kMW0iFRPDY928mfZYPj5AcHDaJKpUqZTXZ758DvsziOzLNOWVhMlArAlTWFg7/PknsHw5k6S//+aHhfWpNJks8Pf3QJ48/FDw8rqOgwdXYeTILqhWzRtVqzKZ8fHh7RMSEpCQkHDvcWJiYlCiRIkMh+QsFn4QHTwIHD8OnDwJnD8PREayYhUdzS/K1IlVrlz8AKldG2jVCmjXTpIoV3X2LLByJbBiBbB5M79ovL2B6tVTEptKlVK+XAICbDuuUnydX7oE/PsvX3979zIJO3WKt8mXD3juOSbprVrJa8xVWZMjaxK+dy9w+jR/FxAAlCnzYBUo9XVwMJPz7EpI4Ofd/dWq1Ana6dM8YQWA0qVTXvuhofwczJ8/538HYVeSMLmC8+eByZOBCRO2I3fuOoiL87mXHHl4cDiscmWgcWNgxYqeaNw4LyZOnHjv/rNmzcKQIUNsbhNgrzlMN28CS5YAa9cCBw4AFy4A8fEpv0+dRPXsCbRoke2HEjo7fhyYP59J0sGDTJCefBJo3Rpo0ACoVo3VTK1Yv0A3b34whrZtgS5d5CzfqJKTgZ07eXJoTZCtyVHu3ClVHGulsXx556gyWizAP/+kxBwRwddo6iTKGvcTT/DiJZNk9CQJk1H9+ScwcSKwaVPq8m8y8uZNRq1afmjcGHjhBeDRR9Peb9iwYfjtt99w6NChez/r0qULbt686RSTvh+WRPn4cD5Uz55A794pFTDhnJKSmJxMmQJs3Ohc1R1rlWvlSr6HvL2Bzp2BgQP5BSWc2+3bwJo1fP5WrQJu3GByVLt2SpJhHYL1MNBab2sSlXqIet8+/nud6f3jpjKfhaaUsvUiNJSUpNTs2Uo1bKiUn59SHIBQKn9+s2rR4qaaNeuEAqC++uortX//fnXu3DmllFLDhw9X3bp1u3ec06dPK39/f/XOO++oY8eOqfDwcOXp6alWr15tcyzR0dEKgIqOjrb7vzM9V64oNXKkUhUqKGUy8d9tMilVtqxS77yj1OXLDglD2OjKFaXGjFGqaFE+Vw0aKDV/vlLx8XpHlr7ISKU+/lipEiUY72OP8b2WmKh3ZCK18+eVCg9XqlUrpXx8+FxVr67Ue+8ptWuXUmaz3hFqw2xWavdupUaNUqpGDf67vb2VatlSqW+/VersWb0jdBuZ5kGSMOnoxg2+SSpVUsrDIyVRKF1aqTffVOrCBaU2bdqkwOpemkuPHj2UUkr16NFDNWnSJM1xN23apGrWrKl8fHxUmTJl1KxZs7IUl6MTptSSkpSaMUOpJ55Qytc3JXEsWFCpTp34wSL0ER3N12tAAC/9+il14IDeUdkuKUmp5cv5hQwwIf/pJ9f9InZ2FotSe/cq9cEHStWqxefEy0upZs2U+uYbpU6f1jtCfZw9y0SpRQsmTgATqVGjlNqzh383oQlJmJzR2rVKhYamVFO8vfn/kyYpdfeu3tHpmzDdb9MmpZ5/Xql8+VKSp0KFlBo9WqmEBL2jcw/x8Up99ZVSBQqw+vnuu0z2jezgQaXatOHrqVYtpVavli8iR4mOVmryZKWqVOHfP29epbp0UWrhQqVu3dI7OucSFaXUokVKde2a8hlYuTK/K6Ki9I7O5UjC5Czu3lVqxAh+6VgrSXXrKrVypfOd4TpTwpTa6dNK9eihVK5cKWejTz+t1NGjekfmupYuVapkSaU8PZXq00epixf1jsi+tm7lkCKg1FNPKXX8uN4Rua6DB5Xq35/VSU9Ppdq3V2rNGhkatVVSklLr1in14ov87AsIUKpvX2NVeZ2cJEx6O3RIqebN+QEB8EXep49zn6E7a8JkZTYrNXOmUmXKpFSdypThUJ6zJZ9Gdf26Up0782/bpo1rJxIWC09cypdnBe2LL5RKTtY7KteQkKDUggWcmwkoFRLC6rCrJd6O9u+/Sn34oXHmERqEJEx6MJs5Bl2yZMoXeoUKSv34o96R2cbZE6bUjh9X6tlnecYFsPrUvTsnJovsWbJEqcKFOQQwf777DFXFxXHuoMnEOXSunCRq7dw5LuQoXJjvy6ZNlfr5Z6km2VtiolKLF3Pel3W6wvDhSp05o3dkhiQJkyOZzTx7sg4ZeXsr1a6d8SYvGilhskpK4hmX9QPa+iEtK+xsFxPDuSSAUm3buu/fbutWpcqVY7Vp0iT3SRjt4Z9/uDjDw0OpPHmUGjxYqSNH9I7KPRw7ptTrrysVGMik/8UXlTpxQu+oDEUSJkf57ju+UAGlgoKUGj+eX+JGMnnyZFW5cmVVoUIFwyVMqW3dmrLqxmRSqkMHpW7f1jsq53b6tFLVqvFLbt48SRLi4pR67TW+hnr2lKGOzFy6pNSAAaz0FivG9gDyntNHbCy/j0qU4FSQvn05hCcyJQmT1pYsUSo4mH9JPz9WOYw+j8aIFab0bN/OSoF1gviAAcZLYh1h40YuRihbVqoB9/vxR7a3qF/ffStuDxMVxaE3f38O4X7+uVJ37ugdlVCKC42+/FKp/Pk56jFsmFI3b+odlVOThEkrf/7JLxjrl/HAga7zZewqCZPVsmVpk9oxY4yf1NpLeDjPQps1c+6FCHratYuTlYsVYx8cwS/jCRNSvoxHjJCWAM4qKorNP/392cLhs88kqc2AJEz2dvx42uGeF190vdKzqyVMVtOmcbjUOmw6bZreEenHYmEXdYDzHlwl2dfKv/+yDUiuXOzZ5K6SkpT6/nulihdnot2/vwz3GMXlyzyx9/Li6rrp0+V9fx9JmOzl7l2lnnkmZULxU0+5boneVRMmpVhZ+vDDlO1nQkLcr3u42czJuIBSEyfqHY1x3Lmj1HPPcduO5cv1jsbxNmxg00RAqY4dlfr7b70jEtlx8mRKy5AKFdhIWSilbMiDDLRtoX5WrAAKFgR+/52
"text/plain": [
"Graphics object consisting of 8 graphics primitives"
]
},
"execution_count": 93,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"draw_circles_from_coords(W) + draw_circles_from_coords(S4 * W)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Calculations to find the curvatures of the dual circles of the octahedral packing\n",
"\n",
"Here the basic idea was to work out the intersection points of the circles and let the computer algebraically find the radii of the circles. The end up coming out, rather anticlimactically, to $(0, 0, \\sqrt 2, \\sqrt 2, \\sqrt 2, \\sqrt 2, 2\\sqrt 2, 2\\sqrt 2)$. It's rather satisfying and a bit surprising that the dual of the root sextuple of the octahedral packing is the root of the cubic packing, since the dual polyhedron of the octahedron is the cube."
]
},
{
"cell_type": "code",
"execution_count": 63,
"metadata": {},
"outputs": [],
"source": [
"def circle_from_points(pta, ptb, ptc):\n",
" a = var('a')\n",
" b = var('b')\n",
" r = var('r')\n",
" x = var('x')\n",
" y = var('y')\n",
" \n",
" circle_func = (x - a)^2 + (y - b)^2 == r^2\n",
" \n",
" eq1 = circle_func.subs(x == pta[0]).subs(y == pta[1])\n",
" eq2 = circle_func.subs(x == ptb[0]).subs(y == ptb[1])\n",
" eq3 = circle_func.subs(x == ptc[0]).subs(y == ptc[1])\n",
" \n",
" res = solve([eq1, eq2, eq3], a, b, r)[1]\n",
" \n",
" return (res[0].rhs(), res[1].rhs(), res[2].rhs())"
]
},
{
"cell_type": "code",
"execution_count": 64,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(1/2*sqrt(2), 1)\n",
"sqrt(2)\n",
"(1/4*sqrt(2), 0)\n",
"2*sqrt(2)\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAGECAYAAAAySIfuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAA9hAAAPYQGoP6dpAABdvElEQVR4nO2dd3gUVRfG300lARIIPRAgEFrohChFEEQRRKpipwiCCCgqKqAi8inSQRGC2MAKKL2qoBBERELoJfSeACEJmxBSd8/3x2GTLCTZNruzuzm/55lnk93ZO2dmZ+57yznnaogIgiAIgmDAQ20DBEEQBOdChEEQBEEwQoRBEARBMEKEQRAEQTBChEEQBEEwQoRBEARBMEKEQRAEQTBChEEQCqBhAjQajUZtWwRBLbys/J5ExQluiVarRWBgILRardqmCIK9MNnokR6DIAiCYIQIgyAIgmCECIMgCIJghAiD4DYsXLgQzZo1Q0BAAAICAtC2bVts3rxZbbMEweUQYRDchho1amDatGmIjY3F3r178dBDD6F37944evSo2qYJgkuhsTLttnglORO3bwNxcUDDhoC/v9rWOBVBQUGYOXMmhg4dWujnWVlZyMrKyvs/NTUVISEh0Gq1CAgIcJSZ7oncl86KeCWVCOLigIgIfhUAADqdDsuWLUN6ejratm1b5H5Tp05FYGBg3hYSEuJAK90cuS9dFhEGwa04fPgwypQpA19fX4wYMQKrV69GeHh4kftPmDABWq02b7t06ZIDrRUE58TaADdBcEoaNGiAAwcOQKvVYsWKFRg0aBCio6OLFAdfX1/4+vo62EpBcG5EGAS3wsfHB2FhYQCAiIgIxMTE4LPPPsOiRYtUtkwQXAcZShLcGr1ebzS5LAiCaaTHILgNEyZMQPfu3VGzZk2kpaXh559/xvbt2/H777+rbZoguBQiDILbcP36dQwcOBAJCQkIDAxEs2bN8Pvvv+ORRx5R2zRBcClEGAS34ZtvvlHbBEFwC2SOQRAEQTBChEEQBEEwQoRBEARBMEKEQRAEQTBChEEQACxYsADh4eGIjIxU2xRBUB0RBkEAMGrUKBw7dgwxMTFqmyIIqiPCIAiCIBghwiAIgiAYIcIgCIIgGCHCIAiCIBghwiAIgiAYIcIgCIIgGCHCIAiCIBghwiAIgiAYIcIgCIIgGCHCIAiCIBghwiAIkFxJglAQEQZBgORKEoSCiDAIgiAIRogwCIIgCEaIMAiCIAhGiDAIgiAIRogwCIIgCEaIMAiCIAhGiDAIgiAIRogwCIIgCEaIMAiCIAhGiDAIAiQlhiAURIRBECApMQShICIMgiAIghEiDIIgCIIRIgyCIAiCESIMgiAIghEiDIIgCIIRIgyCIAiCESIMgiAIghEiDIIgCIIRIgyCIAiCESIMgiAIghEiDIIAyZUkCAURYRAESK4kQSiICIMgCIJghAiDIAiCYIQIgyAIgmCECIMgCIJghAiDIAiCYIQIgyAIgmCECIPgNkydOhWRkZEoW7YsKleujD59+uDEiRNqmyUILocIg+A2REdHY9SoUdi9eze2bNmCnJwcdO3aFenp6WqbJgguhZfaBgiCUvz2229G/y9ZsgSVK1dGbGwsOnbsqJJVguB6iDAIbotWqwUABAUFFblPVlYWsrKy8v5PTU21u12C4OzIUJLgluj1erz++uto3749mjRpUuR+U6dORWBgYN4WEhLiQCsFwTkRYRDcklGjRuHIkSNYtmxZsftNmDABWq02b7t06ZKDLBQE50WGkgS3Y/To0diwYQN27NiBGjVqFLuvr68vfH19HWSZILgGIgyC20BEePXVV7F69Wps374doaGhapskCC6JCIPgNowaNQo///wz1q5di7Jly+Lq1asAgMDAQPj5+alsnSC4DjLHILgNCxcuhFarRadOnVCtWrW8bfny5WqbJgguhfQYBLeBiNQ2QRDcAukxCIIgCEaIMAiCIAhGiDAIgiAIRogwCIIgCEaIMAiCIAhGiDAIAoAFCxYgPDwckZGRapsiCKojwiAI4OC4Y8eOISYmRm1TBEF1RBgEQRAEI0QYBEEQBCMk8llQntxc4No1ID4eSEgAUlKAnByACPDyAvz8gKpVgeBgoFo1oGxZtS12XoiApCS+jvHxwPXrfC11OsDTE/DxASpVyr+WFSoAGo3aVgsujgiDYBtpacCePUBsbP529ixXaOYSGAi0bAlERPB2331A3br2s9lZyc0FDhwA9u7Nv5ZHjwLZ2eaX4esLhIfnX8vWrYEWLViQBcFM5G4RLOfyZWD9emDdOuCvv7jiKlOGK/eePbliMrRgg4OBoCCumDQaQK8H0tO5BWzYzp0D9u0DVq4EZs/mY9SvD/TqxVvbtu5bsaWlAb//ztdz40buHXh6Ak2acMU+aBAQEpJ/LStX5srfw4OvZVYW9yIMvbNLl4BDh1isFy/mnkWFCkCPHnwtu3aVHppgEjd92gTFyc4GVq8GoqKAHTu4on7wQWDmTOCRR4AGDbiyMoWnJxAQwFuDBvd+npQE/PMPV5Q//ADMmsWV4UsvAS+/DNSsqfy5ORoiYPduvpa//sqVe9OmfH6PPQa0asXDbabw8OD9atXi7W4yMlhwN21iEf/+exaV/v2BkSOBNm1k2EkoHCKyZhOcidhYIoBflSYxkej994mqVOFjPPgg0Y8/EqWkKH+su9HpiHbvJho9migggMjDg6hXL6Lt2+12SK1WSwBIq9UqX3h2NtHXXxO1aMHXsm5dohkziM6eVf5YhXH2LNH06UR16vDxW7Rge7Kz7XM8e96Xgi2YrONFGNwBezyAaWlEkycTlS1LVKYM0ahRREeOKFe+NfYsWkTUtCmf66OPEu3bp/hh7CIMOh3RsmVEYWFEGg1Rz55Emzfz+2qg0/Hxe/Zke8LC2D6l7RFhcFZM1vHirioYo9cDCxfy5O+UKcCwYTwHMH8+0LixenaVKQMMHw4cPAisWAGcP89DLs8+y3MezsqOHUBkJPDMMzx0duAAD+t062be0Js98PDg469bB+zfz/M5zzzDdu7YoY5NglMhwiDkc/o00KkTjz937w6cOsWTwRUrqm1ZPhoN8MQTwJEjwFdfAdHRLFjffmuZJ5S9uXULePVVnofx9uYKd8MGoFkztS0zpnlznvSOjmY7H3wQeO01dhAQSiwiDAL3EubN40rr8mX2NFqyxLkner28eEL66FGgb19g6FCeuLWy96BorqToaL6W33wDfPopsGsX0KGD7eXak44d2c5PPwW+/prtl95DiUWEoaRz6xbw5JPAmDHAkCHs6ti5s9pWmU/58ixiGzaw7a1aAX//bXExiuRKIuIe1kMPAdWr87DXmDHqDRlZiocH23vwILvGdu4MzJnjXD0xwSG4yB0r2IXz54F27YAtW4A1a3geoUwZta2yjh49uEJr3Bjo0oWHmRxJZiYweDDw1lvA228D27cD9eo51galqFeP7X/7bWDsWODFF/n8hBKDCENJZdcunmxMT2ef+t691bbIdipWBP74gyfMhw/n1q9eb//j3rjBczO//AL89BMwbRrHa7gynp58Hj/+CCxfzr2HGzfUtkpwECIMJZHt2zkCNjycI2TV9DZSGm9vYMEC9qz6/HOee9Dp7He8q1dZFM6d4yGs556z37HU4Pnnea7h7Fk+z2vX1LZIcAAiDCWNHTt4krZdO2DzZk6X4I6MGMGt9x9+4LkTe/Qcrl/n+YSUFJ5wbt1a+WM4A5GRfH4pKdxzuH5dbYsEOyPCUJLYu5fH4tu3Zx92f3+1LbIvzz7LwvDDD+w6quQkamoqpwJJSQG2bQMaNlSubGekYUM+z+RkPu/UVLUtEuyICENJIT6e5xGaNOGJ5lKl1LbIMTz7LPDll5yXKCpKmTJ1Oh4yunAB2LqVA8RKAvXrA3/+yef9/PP2HaITVEWEoSSQmcm+/hoNJ8IrXVptixzLSy8Br7/Ok9F//WV7ee+/z8Nwy5a51/yMOTRuDCxdykFxEyeqbY1gJ0QYSgIvv8w+/mvW8AI5JZGZM3k+oH9/nki1lqVL2VtnxgxOK1ES6d6dz3/
"text/plain": [
"Graphics object consisting of 14 graphics primitives"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"c1 = line([(-3, 1), (3, 1)])\n",
"c2 = line([(-3, -1), (3, -1)])\n",
"c3 = circle((-sqrt(2), 0), 1)\n",
"c4 = circle((sqrt(2), 0), 1)\n",
"c5 = circle((0, 1/2), 1/2)\n",
"c6 = circle((0, -1/2), 1/2)\n",
"\n",
"b1 = line([(-sqrt(2), -3), (-sqrt(2), 3)], rgbcolor=(1, 0, 0))\n",
"b2 = line([(sqrt(2), -3), (sqrt(2), 3)], rgbcolor=(1, 0, 0))\n",
"\n",
"x, y, r = circle_from_points((sqrt(2) / 3, 1 / 3), (0, 1), (sqrt(2), 1))\n",
"print('({}, {})'.format(x, y))\n",
"print(1 / r)\n",
"\n",
"b3 = circle(( x, y), r, rgbcolor=(1, 0, 0))\n",
"b4 = circle((-x, y), r, rgbcolor=(1, 0, 0))\n",
"b5 = circle(( x, -y), r, rgbcolor=(1, 0, 0))\n",
"b6 = circle((-x, -y), r, rgbcolor=(1, 0, 0))\n",
"\n",
"x, y, r = circle_from_points((sqrt(2) / 3, 1 / 3), (0, 0), (sqrt(2) / 3, -1 / 3))\n",
"print('({}, {})'.format(x, y))\n",
"print(1 / r)\n",
"\n",
"b7 = circle(( x, y), r, rgbcolor=(1, 0, 0))\n",
"b8 = circle((-x, y), r, rgbcolor=(1, 0, 0))\n",
"\n",
"show(c1 + c2 + c3 + c4 + c5 + c6 + b1 + b2 + b3 + b4 + b5 + b6 + b7 + b8)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Futile attempt to derive the quadratic form for the cubic packing\n",
"\n",
"We attempted to derive the quadratic form for the cubic packing based on the data collected above, but the results don't match the results from the Stange at all and are significantly uglier. There's probably a mistake in here somewhere, but I'm not sure where."
]
},
{
"cell_type": "code",
"execution_count": 65,
"metadata": {},
"outputs": [],
"source": [
"def abbc_coords(b, h1, h2):\n",
" return [(h1^2 + h2^2 - 1) / b, b, h1, h2]"
]
},
{
"cell_type": "code",
"execution_count": 66,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[ 1 -3 -3 -5/4 -1 -1 -5 3/4]\n",
"[ -3 1 -3 1/2 -1 -1 -1 -3/2]\n",
"[ -3 -3 1 -5/4 -5 -1 -1 3/4]\n",
"[-5/4 1/2 -5/4 1 3/4 -3/2 3/4 -1]\n",
"[ -1 -1 -5 3/4 1 -3 -3 -5/4]\n",
"[ -1 -1 -1 -3/2 -3 1 -3 1/2]\n",
"[ -5 -1 -1 3/4 -3 -3 1 -5/4]\n",
"[ 3/4 -3/2 3/4 -1 -5/4 1/2 -5/4 1]"
]
},
"execution_count": 66,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Wc = matrix([\n",
" abbc_coords(sqrt(2), -1, sqrt(2)),\n",
" abbc_coords(2*sqrt(2), 1, 0),\n",
" abbc_coords(sqrt(2), -1, -sqrt(2)),\n",
" [sqrt(2) / 4, 0, 1, 0],\n",
" abbc_coords(sqrt(2), 1, sqrt(2)),\n",
" abbc_coords(2*sqrt(2), -1, 0),\n",
" abbc_coords(sqrt(2), 1, -sqrt(2)),\n",
" [sqrt(2) / 4, 0, -1, 0]\n",
"\n",
"])\n",
"\n",
"Wc * P * Wc.transpose()"
]
},
{
"cell_type": "code",
"execution_count": 67,
"metadata": {},
"outputs": [],
"source": [
"W = matrix([\n",
" abbc_coords(sqrt(2), -1, sqrt(2)),\n",
" abbc_coords(2*sqrt(2), 1, 0),\n",
" abbc_coords(sqrt(2), -1, -sqrt(2)),\n",
" [sqrt(2) / 4, 0, 1, 0]\n",
"])"
]
},
{
"cell_type": "code",
"execution_count": 68,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/plain": [
"[ 1 -3 -3 -5/4]\n",
"[ -3 1 -3 1/2]\n",
"[ -3 -3 1 -5/4]\n",
"[-5/4 1/2 -5/4 1]"
]
},
"execution_count": 68,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"m = W * P * W.transpose()\n",
"m"
]
},
{
"cell_type": "code",
"execution_count": 69,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"97*b1^2 - 304*b1*b2 + 328*b2^2 - 290*b1*b3 - 304*b2*b3 + 97*b3^2 + 32*b1*b_avg - 1088*b2*b_avg + 32*b3*b_avg + 1280*b_avg^2"
]
},
"execution_count": 69,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"D = W2.transpose() * m.inverse() * W2\n",
"968*factor(simplify(D[1][1]))"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "SageMath 9.2",
"language": "sage",
"name": "sagemath"
},
"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.2"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Reference in a new issue Copy permalink