Hexagonal number spirals

Having recently discovered that the prevailing number spirals of our age are not nearly hexagonal enough for my needs, I invented new a circular spiral based on centered hexagonal numbers, which I have named the GRAHAM SPIRAL. The first rotation contains 6 dots, representing the first 6 natural numbers, and each subsequent rotation n contains n*6 dots, in the fashion of hexagonal tessellation. In the examples above, variable dot sizes represent the number of unique prime factors in a given number, while the uniform dots represent primes themselves, but of course any arithmetic attributes can be plotted in the same fashion. The resulting patterns clearly demonstrate the deeply hexagonal structure of the natural numbers, &c.

(Style note: "GRAHAM SPIRAL," in sentence form, should always be written in ALL CAPS.)

I came upon for the idea for this spiral while investigating the Sacks spiral, which is based on one perfect square for each full rotation. That is, its rotations proceed 4, 9, 16... instead of 6, 12, 18... While certainly better (and indeed more hexagonal) than the square Ulam spiral, it still suffers from some obvious deficiencies. Circles don't pack around other circles in radiating squares, they pack in hexagons. Thus, a hexagonal spiral ought to be considered the "natural" way to plot a circular number spiral on the Euclidean plane. I mean, I'm not saying the Sacks spiral isn't interesting, or wouldn't be appropriate in another topological context, I'm just saying it doesn't exactly accentuate the organizational advantages of the Euclidean plane the way a hexagonal spiral does. I think there is a natural intuitiveness and ergonomicity to the GRAHAM SPIRAL that gives it distinct advantages over other, lesser spirals.

It doesn't particularly matter where you start, as it will eventually even out into roughly the same pattern, though you can twist it either way by adjusting the starting angle. In these examples, the 1-dot is plotted at approximately sixty degrees in the top right (positive) quadrant, spiraling counterclockwise, which puts 6, 60, and other highly composite numbers on the straight line going along the positive X axis from the origin. The 1, 7, 19... line (i.e., the proper centered hexagonal numbers) have more primes, so if you're into primes or something, it might make sense to center that line instead. This could also be accomplished by putting the 1 in the center, instead of the theoretical zero.

Anyway, The Python script—which I adapted from the Sacks spiral script at http://www.dcs.gla.ac.uk/~jhw/spirals/—is below, so the interested reader may recreate their own GRAHAM SPIRALS as they see fit. I edited this in a disjointed, nonlinear fashion over the course of several days, so it could probably stand to be cleaned up a bit. It requires the Elementary Number Theory (ent) package, which I can't find at the moment, but which I assure you is somewhere on these very internets at this time. Thank you.

# GRAHAM SPIRAL: The official script
from pyx import *
from ent import *
from math import *

ca = canvas.canvas()
tau = 2 * pi
off = 0
n = input("Rings:")
n = (n+1)*(n+1)*(n+1) - n*n*n - 1

# Use this to add an arbitary numerical offset
#off = input("Offset:")

# What "ring" any given natural number falls into
def hexring(n):
    val = 0
    while n > 0:
        val += 1
        n = n - 6 * val
    return val

# The nth hex number - 1, based on being the gnomon of the nth cube
def hexnum(n):
    return (n+1)*(n+1)*(n+1) - n*n*n - 1

for j in range(n):
    i = j + 1
    ring = hexring(i) # What ring is i in
    perim = ring * 6 # How many numbers are in the perimeter of that ring
    ringcount = hexnum(ring) - i # Where is i in that ring
    r = ring + 1 - (float(1) / float(perim)) * ringcount
    theta = (tau / perim) * ringcount
    x = cos(theta) * r
    y = -sin(theta) * r
    factors = factor(i+off)
    
    # Uncomment this to produce dotted primes
    #numfac = 0
    #for f in factors:
    #    numfac += f[1]
    #if numfac == 1: radius = 0.4
    #else: radius = 0.1
    
    # Set dot size by number of unique prime factors
    npf = len(factors)
    radius = pow(1.414, npf) / 12

    ca.fill(path.circle(x, y, radius))

d = document.document(pages = [document.page(ca, paperformat=document.paperformat.A4, fittosize=1)])
d.writePSfile("grahamspiral.ps")