Snub-dodecahedron-python
Jump to navigation
Jump to search
Numerical Solution to Snub Dodecahedron
#========================================
# Numerical Solution to Snub Dodecahedron
#========================================
# mpmath is included in Anaconda and Sympy
# Set accuracy to 50 decimal places
import mpmath as mp
mp.mp.dps = 50
class SnubNumerical(object):
def __init__(self, log):
self.log= log
class Point(object):
def __init__(self, parent, x, y, z):
self.parent= parent
self.x= x
self.y= y
self.z= z
def make_copy(self):
self.mcopy= self.parent.Point(self.parent,self.x,self.y,self.z)
return self.mcopy
def make_vector_to(self, endpoint):
self.vector = self.parent.Point(self.parent,
endpoint.x-self.x,
endpoint.y-self.y,
endpoint.z-self.z)
return self.vector
def distance_to(self, endpoint):
self.distance = mp.sqrt(
(endpoint.x-self.x)**2 +
(endpoint.y-self.y)**2 +
(endpoint.z-self.z)**2 )
return self.distance
def add_vector(self, point):
self.x += point[0]
self.y += point[1]
self.z += point[2]
return
def scale(self, ratio):
self.ratio = ratio
self.scaled_xyz = [
self.x * self.ratio,
self.y * self.ratio,
self.z * self.ratio]
return self.scaled_xyz
class Triangle(object):
def __init__(self, parent, pointA, pointB, pointC):
self.parent= parent
self.point_A= parent.Point(self.parent, pointA[0], pointA[1], pointA[2])
self.point_B= parent.Point(self.parent, pointB[0], pointB[1], pointB[2])
self.point_C= parent.Point(self.parent, pointC[0], pointC[1], pointC[2])
self.center = parent.Point(self.parent,
(self.point_A.x + self.point_B.x + self.point_C.x) /3,
(self.point_A.y + self.point_B.y + self.point_C.y) /3,
(self.point_A.z + self.point_B.z + self.point_C.z) /3)
self.ratio_1= mp.mp.mpf('.1')
self.ratio_2= mp.mp.mpf('.1')
self.vector_to_g_1 = self.center.make_vector_to( self.point_B)
self.vector_to_g_2 = self.point_B.make_vector_to(self.point_C)
self.vector_to_j_1 = self.center.make_vector_to( self.point_C)
self.vector_to_j_2 = self.point_C.make_vector_to(self.point_A)
def set_ratio_1(self, ratio_1):
self.ratio_1= ratio_1
return
def set_ratio_2(self, ratio_2):
self.ratio_2= ratio_2
return
def calculate(self):
self.point_g = self.center.make_copy()
self.point_g.add_vector(self.vector_to_g_1.scale(self.ratio_1))
self.point_g.add_vector(self.vector_to_g_2.scale(self.ratio_2))
self.point_j = self.center.make_copy()
self.point_j.add_vector(self.vector_to_j_1.scale(self.ratio_1))
self.point_j.add_vector(self.vector_to_j_2.scale(self.ratio_2))
# Ratio from center of eq. triangle to vertex by edge length is sqrt(3)
self.distance_g_j = self.point_g.distance_to(self.center) * mp.sqrt(3)
return
def info(self):
digits = 5
D1 = mp.nstr(self.point_A.distance_to(self.point_B), digits)
D2 = mp.nstr(self.point_B.distance_to(self.point_C), digits)
D3 = mp.nstr(self.point_C.distance_to(self.point_A), digits)
self.parent.log.info('Base: %s\n %s\n %s'%(D1,D2,D3))
return
def run_numerical_solution(self, use_snub_cube=False):
self.use_snub_cube = use_snub_cube
mp.mp.dps = 50
self.verbose = True
if self.use_snub_cube:
self.make_icosa()
self.triangle_1 = self.Triangle(self, self.icosa[1], self.icosa[0], self.icosa[2])
self.triangle_2 = self.Triangle(self, self.icosa[3], self.icosa[2], self.icosa[0])
else:
self.make_icosa()
self.triangle_1 = self.Triangle(self, self.icosa[1], self.icosa[0], self.icosa[2])
self.triangle_2 = self.Triangle(self, self.icosa[3], self.icosa[2], self.icosa[0])
self.finder_counter = 0
tolerance = mp.mpf('1.0e-33')
mp.findroot(self.finder_f, mp.mpc('0.1','0.1'), tol=tolerance,
solver='halley', maxsteps=2000, verbose=False)
if self.verbose:
if self.use_snub_cube:
self.log.info("Numerical solution to Snub Cube")
else:
self.log.info("Numerical solution to Snub Dodecahedron")
self.log.info("Findroot (halley) ran %d iterations, tolerance %s"%(
self.finder_counter, mp.nstr(tolerance, 7)))
digits = 31
self.log.info(' v v v')
self.log.info('Segment g_j: = %s'%(mp.nstr(self.triangle_1.distance_g_j, digits)))
self.log.info('Segment g_gprime = %s'%(mp.nstr(self.distance_g_gprime, digits)))
self.log.info('Segment j_gprime = %s'%(mp.nstr(self.distance_j_gprime, digits)))
if self.use_snub_cube:
self.log.info('snub_cube not done')
else:
phi = (mp.sqrt(5) + 1) / 2 # Golden Section
sqr1 = mp.sqrt( phi**2 / 4 - mp.mpf(8) / 27 )
xi = mp.exp( mp.log( phi / 2 + sqr1) / 3) + mp.exp( mp.log( phi / 2 - sqr1) / 3)
D = phi /(mp.sqrt( phi**2 + 3 * xi * (phi + xi)))
self.log.info('Calculated D = %s'%(mp.nstr(D, digits)))
return
def finder_f(self, ratio):
'Called by mp.findroot'
self.finder_counter += 1
self.triangle_1.set_ratio_1(ratio.real)
self.triangle_2.set_ratio_1(ratio.real)
self.triangle_1.set_ratio_2(ratio.imag)
self.triangle_2.set_ratio_2(ratio.imag)
self.triangle_1.calculate()
self.triangle_2.calculate()
self.distance_j_gprime = self.triangle_1.point_j.distance_to(self.triangle_2.point_g)
self.distance_g_gprime = self.triangle_1.point_g.distance_to(self.triangle_2.point_g)
self.delta_1 = self.triangle_1.distance_g_j - self.distance_j_gprime
self.delta_2 = self.triangle_1.distance_g_j - self.distance_g_gprime
self.delta_distance = mp.mpc(self.delta_1, self.delta_2)
return(self.delta_distance)
def make_icosa(self, verbose=0, unit_edge_length=1 ):
self.icosa_faces = [
( 0, 1, 2) ,( 0, 2, 3) ,( 0, 3, 4) ,( 0, 4, 5) ,( 0, 5, 1),
(11, 6, 7) ,(11, 7, 8) ,(11, 8, 9) ,(11, 9, 10) ,(11, 10, 6),
( 1, 2, 6) ,( 2, 3, 7) ,( 3, 4, 8) ,( 4, 5, 9) ,( 5, 1, 10),
( 6, 7, 2) ,( 7, 8, 3) ,( 8, 9, 4) ,( 9, 10, 5) ,(10, 6, 1)]
c63 = 1 / mp.sqrt(5) # Cos(a) a = arctan(2) = 63.434... degrees
s63 = 2 / mp.sqrt(5) # Sin(a) a = arctan(2) = 63.434... degrees
c72 = mp.sqrt((3-mp.sqrt(5))/8) # Cos(72)
s72 = mp.sqrt((5+mp.sqrt(5))/8) # Sin(72)
c36 = mp.sqrt((3+mp.sqrt(5))/8) # Cos(36)
s36 = mp.sqrt((5-mp.sqrt(5))/8) # Sin(36)
self.icosa = [
[ 0, 0, 1],
[ s63 , 0, c63],
[ s63*c72, s63*s72, c63],
[-s63*c36, s63*s36, c63],
[-s63*c36,-s63*s36, c63],
[ s63*c72,-s63*s72, c63],
[ s63*c36, s63*s36,-c63],
[-s63*c72, s63*s72,-c63],
[-s63 , 0,-c63],
[-s63*c72,-s63*s72,-c63],
[ s63*c36,-s63*s36,-c63],
[ 0, 0, -1]]
if unit_edge_length:
self.point_A= self.Point(self, self.icosa[0][0],self.icosa[0][1],self.icosa[0][2])
self.point_B= self.Point(self, self.icosa[1][0],self.icosa[1][1],self.icosa[1][2])
adjust = 1 / self.point_A.distance_to(self.point_B)
for i, points in enumerate(self.icosa):
for j, xyz in enumerate(points):
self.icosa[i][j] *= adjust
if verbose:
digits = 3
for i, xyz in enumerate(self.icosa):
x = mp.nstr(xyz[0], digits)
y = mp.nstr(xyz[1], digits)
z = mp.nstr(xyz[2], digits)
print('V %d (%s, %s %s)'%(i+1, x,y,z))
return
if __name__ == '__main__':
snub_numerical = SnubNumerical(MyLog())
snub_numerical.run_numerical_solution()
notes = '''
Output:
Numerical solution to Snub Dodecahedron
Findroot (halley) ran 492 iterations, tolerance 1.0e-33
v v v
Segment g_j: = 0.3638558935244065463166939410955
Segment g_gprime = 0.3638558935244065463166939410955
Segment j_gprime = 0.3638558935244065463166939410955
Calculated D = 0.3638558935244065463166939410955
'''
#===============================================
# End of Numerical Solution to Snub Dodecahedron
#===============================================