-
Notifications
You must be signed in to change notification settings - Fork 3
/
ammonite.py
191 lines (153 loc) · 6.73 KB
/
ammonite.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
# Generate a snail (ammonite) shell from triangles
# This method is based in the one decribed in Spirals
# by Tomoko Fuse for Ammonites.
# The flat method is controlled by a cental angle (ca) and
# an angle of spirality (sa). With a starting length (u), and a
# number of iterations (N), this is enough to create
# the ammonite pattern.
# Start with an isoceles triangle where the central angle (ca)
# is the one that is the odd one out. The starting length is the
# leg opposite ca and is 2u. Drop a bisector from ca to get a pair
# of right triangles. This is the tip of the ammonite.
# Next, add a triangle from the center line (bisector) such
# that it has an angle sa from u. Extend this and extend the legs
# of the isosceles triangle until they meet.
# Add a line from this point to the center line (perpindicular
# to the centerline), this becomes u_prime
#
# The resulting shape is a right trapezoid with a diagonal.
# Mirror this on the other side of the centerline.
# Set u=u_prime and repeat N times.
#
# While the diagonal cuts the trapezoid into a right
# triangle and a scalene trianle, sadly, we know information
# for the scalene trianle but need information for the right
# triangle which means we need
#
# The Law of Cosines: Pythagorean theorem for non-right triangles
# a*a = b*b + c*c - 2ac cos(angle_a)
# The Law of Sines (this is far more useful for this)
# a/sin(angle_a) = b/sin(angle_b) = c/sin(angle_c)
# Much algebra will be done. I will not show my work because
# I mostly used a website to get the equations
# https://www.mathsisfun.com/algebra/trig-solving-triangles.html
# And anothe website to check my work:
# https://www.triangle-calculator.com/
#
# A straight line approach leads to flat origami. I modified
# the algorithm to use bezier curves to get a self-shaping
# curved shell.
from math import *
import matplotlib.pyplot as plt
import numpy as np
from matplotlib import collections as mc
from matplotlib.patches import Ellipse, Wedge, Polygon
import itertools
class point:
def __init__(self, x, y):
self.x = x
self.y = y
def lengthTo(self, p):
x2 = (self.x-p.x)**2
y2 = (self.y-p.y)**2
len = sqrt(x2 + y2)
return len
def pts(self):
return [self.x, self.y]
def negative(self):
return point(self.x, -self.y)
def pointFrom(self, length, angle):
x = self.x + length*sin(radians(angle))
y = self.y + length*cos(radians(angle))
return point(x, y)
class side:
def __init__(self, pointA, pointB):
self.pointA = pointA
self.pointB = pointB
def length(self):
return self.pointA.lengthTo(self.pointB)
def add_plot(ax, pt1, pt2, ptb, color):
bezier(ax, pt1, pt2, ptb, color)
def bezier(ax, pt1, pt2, ptb, color):
npoints = 2
numbers = [i for i in range(npoints)]
bezier_path = np.arange(0.0, 1.01, 0.01)
x1y1 = x1, y1 = pt1.x, pt1.y
x2y2 = x2, y2 = pt2.x, pt2.y
xbyb = xb, yb = ptb.x, ptb.y
# Compute and store the Bezier curve points
x = (1 - bezier_path)** 2 * x1 + 2 * (1 - bezier_path) * bezier_path * xb + bezier_path** 2 * x2
y = (1 - bezier_path)** 2 * y1 + 2 * (1 - bezier_path) * bezier_path * yb + bezier_path** 2 * y2
ax.plot(x, y, color)
# We are given two angles of a triangle and one side, which is not the side adjacent to the two given angles. Return all three point coordinates.
# https://www.mathsisfun.com/algebra/trig-solving-triangles.html
# https://www.triangle-calculator.com/
def triangle_find_vertices_aas(angleA, angleC, sideAB):
lengthBC= sideAB.length()*sin(radians(angleA))/sin(radians(angleC))
angleB = 180 - abs(angleA) - abs(angleC)
#print(angleA, angleB, angleC)
lengthAC = sideAB.length()*sin(radians(angleB))/sin(radians(angleC))
pointC = sideAB.pointA.pointFrom(lengthAC, angleA)
#print('AC len {} BC len {} AB len {} BC {} AC {}'.format(lengthAC, lengthBC, sideAB.length(), pointC.lengthTo(sideAB.pointB), pointC.lengthTo(sideAB.pointA)))
return [sideAB.pointA, sideAB.pointB, pointC]
def ptb_straightline(pt1, pt2):
return pt1
def make_plot(ax, ca=20, sa=15, u=1.0, N=15, curve_fun=ptb_straightline, distmult=1):
angleC = 90.0 - ca/2.0 - sa
angleA = sa
pointA = point(0,0)
pointB = pointA.pointFrom(u, 0)
sideAB = side(pointA, pointB)
pointA, pointB, pointC = triangle_find_vertices_aas(angleA, angleC, sideAB)
outerPolyVertices = [pointC.pts(), pointA.pts(), pointC.negative().pts()]
lengthAC = pointA.lengthTo(pointC)
lengthANewA = lengthAC * sin(radians(angleA))
pointA = pointA.pointFrom(lengthANewA*distmult, 90)
pointB = pointC.pointFrom(lengthANewA*(distmult-1.0), 90)
for i in range(N):
sideAB = side(pointA, pointB)
pointA, pointB, pointC = triangle_find_vertices_aas(angleA, angleC, sideAB)
add_plot(ax, pointA, pointB, color='r',
ptb=curve_fun(pointA, pointB))
add_plot(ax, pointA, pointB.negative(), color='r',
ptb=curve_fun(pointA, pointB.negative()))
add_plot(ax, pointA, pointC, color='b',
ptb=curve_fun(pointA, pointC))
add_plot(ax, pointA, pointC.negative(), color='b',
ptb=curve_fun(pointA, pointC.negative()))
lengthAB = sideAB.length()
lengthAC = pointA.lengthTo(pointC)
newLengthAB = lengthAC * sin(radians(90-angleA))
lengthANewA = lengthAC * sin(radians(angleA))
pointA = pointA.pointFrom(lengthANewA*distmult, 90)
pointB = pointC.pointFrom(lengthANewA*(distmult-1.0), 90)
# finish making outer triangle for cutting
sideAB = side(pointA, pointB)
pointA, pointB, pointC = triangle_find_vertices_aas(angleA, angleC, sideAB)
outerPolyVertices.extend([pointC.negative().pts(), pointC.pts()])
poly = Polygon(outerPolyVertices, facecolor='1.0', edgecolor='k')
ax.add_patch(poly)
plt.axis('off')
ax.set_aspect(1), ax.autoscale()
# example values
ca = 30.0 # central angle, book uses 30 to 45 degrees
sa = 15.0 # spirality angle book uses 12 to 15 degrees
u = 1.0 # initial length
N = 20 # number of iterations
show_plot = True
figure, ax = plt.subplots()
name = 'fg24_ca30_sa15_N20'
make_plot(ax, ca=30, sa=15, N=20, curve_fun=ptb_straightline)
plt.savefig(name + ".svg")
if show_plot: plt.title(name), plt.show()
figure, ax = plt.subplots()
name = 'fg28_ca30_sa10_N20'
make_plot(ax, ca=30, sa=10, N=20, curve_fun=ptb_straightline)
plt.savefig(name + ".svg")
if show_plot: plt.title(name), plt.show()
figure, ax = plt.subplots()
name = 'fg28_ca45_sa15_N12_dist2'
make_plot(ax, ca=30, sa=15, N=12, curve_fun=ptb_straightline, distmult=2)
plt.savefig(name + ".svg")
if show_plot: plt.title(name), plt.show()
print(".")