gosper.py

from __future__ import division
import numpy as np

import matplotlib.pyplot as plt

from copy import copy

"""
This script draws the Gosper island, a hexagon-based fractal that tiles the plane.
https://en.wikipedia.org/wiki/Gosper_curve#Properties
"""

n = 5

def hexwalk(n):
   if n == 1:
       return 'FR'*6
   return hexwalk(n-1).replace('F', 'FRFLF')

radius = 100.
mass = 1.

step_size = radius / 7**(.5*(n-1))

cur_angle = 0.0
cur_pos = np.array([0.0, 0.0])
my_positions = []
positions = []

for char in hexwalk(n):
    if char == 'F':
        cur_pos += np.array([step_size*np.cos(cur_angle), step_size*np.sin(cur_angle)])
    elif char == 'R':
        cur_angle -= np.radians(60.0)
    elif char == 'L':
        cur_angle += np.radians(60.0)
    my_positions.append(copy(cur_pos))

my_positions = np.vstack(my_positions)

plt.plot(my_positions[:, 0], my_positions[:, 1])
plt.show()

unique_positions = []
cur_row = None
for row in my_positions:
    if cur_row is None:
        cur_row = row
        unique_positions.append(row)
    else:
        if (row == cur_row).all():
            continue
        else:
            unique_positions.append(row)
            cur_row = row

unique_positions = np.vstack(unique_positions)
unique_positions -= np.sum(unique_positions, axis=0)/unique_positions.shape[0]

#angles = np.linspace(-np.pi, np.pi, 200)
#unique_positions = np.vstack((radius*np.cos(angles), radius*np.sin(angles))).T

plt.plot(unique_positions[:, 0], unique_positions[:, 1])
plt.show()

def rot(theta):
    return np.array([[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]])

def inertia_triangle(p1, p2):
    """Polar inertia around the z axis for a triangle formed by p1, p2 and the
    origin """
    a = np.linalg.norm(p2 - p1)
    b = np.linalg.norm(p1)
    c = np.linalg.norm(p2)

    return -1/48*np.sqrt(-(a-b-c)*(a+b-c)*(a-b+c)*(a+b+c))*(a**2-3*(b**2+c**2))

def area_triangle(p1, p2):
    a = np.linalg.norm(p2 - p1)
    b = np.linalg.norm(p1)
    c = np.linalg.norm(p2)

    return 1/4*np.sqrt((a+b+c)*(b+c-a)*(c+a-b)*(a+b-c))

inertia = 0
area = 0

for i in range(unique_positions.shape[0]):
    #tri = plt.Polygon([np.array([0, 0]), unique_positions[i-1, :], unique_positions[i, :]])
    #fig = plt.figure()
    #ax = fig.add_subplot(111)
    #ax.plot(unique_positions[:, 0], unique_positions[:, 1])
    #ax.add_patch(tri)
    #plt.show()
    inertia += inertia_triangle(unique_positions[i-1, :], unique_positions[i, :])
    area += area_triangle(unique_positions[i-1, :], unique_positions[i, :])

true_m = mass*radius**2.*3./4
print("Radius : {}".format(radius))
print("Total inertia: {}".format(inertia))
print("Inertia/area : {}".format(inertia/area))
print("Error: {}".format(abs((inertia/area - true_m)/true_m)))