```"""
by K. Urner, 4D Solutions
Oct 4, 2000:    fixed jitterbug.drawrect to allowing edge filtering (a feature
dropped from some earlier version -- e.g. see April 13)
May 9, 2000:    moved to Qvectors for default tetrahedron
added rotation method, split out jitterbug-related methods
and moved to multiple inheritance for the relevant shapes,
added Cube subclass
May 7, 2000:    make volume methods work with quadrays
May 4, 2000:    added volume3 method
April 30, 2000: added new tetrahedron volume method
April 22, 2000: added _div_ method for shape/scalar operation
April 13, 2000: added self.control to Shape object
(for drawfilter to work right)
April 3, 2000:
Added scale function to all shapes, enhanced volume()
capabilities accordingly, gave Tetra a vertdict so it
would work in print utilities
Mar 28, 2000:
added print utilities
last modified Mar 23, 2000
minor tweek to Octa to better conform with essay
March 5:
streamlined Tetra class (no arguments required)
used Tetra volume method to derive Icosahedron's
"""

from povray import Povray
from coords import Vector, Qvector, length
import functions
import quat
import rand
import copy

phi = (1 + 5**0.5)/2.0
root2 = 2.0**0.5
syn3 = (9.0/8.0)**0.5

a = Qvector((1,0,0,0))
b = Qvector((0,1,0,0))
c = Qvector((0,0,1,0))
d = Qvector((0,0,0,1))

class Shape:

control = 1.0  # edge length

vertdict = {}
edgepairs = []

def getverts(self,keys):
# utility to retrieve a list of vertices
return map(lambda x,dict=self.vertdict: dict[x],keys)

def __mul__(self,scalefactor):
for i in self.vertdict.keys():
self.vertdict[i] = self.vertdict[i] * scalefactor
self.control = self.control * scalefactor

__rmul__ = __mul__

def __div__(self,scalefactor):
self.__mul__(1.0/scalefactor)

def rotate(self, axis, alpha):
# axis is 'X','Y' or 'Z'
# alpha in degrees
for i in self.vertdict.keys():
self.vertdict[i]=quat.rotate(self.vertdict[i],axis,alpha)

def drawself(self,myfile):
# draw the edges defined in the edgepairs list
for i in self.edgepairs:
myfile.edge(self.vertdict[i[0]],self.vertdict[i[1]])

class Jitterbug:

# Edges of 8 triangles defined by the 12 vertices generated
# as the corners of (a) root2 x root2 squares (cuboctahedron)
# or (b) 1 x phi rectangles (icosahedron)
edgepairs = [['Axy','Axz'],['Axz','Ayz'],['Ayz','Axy'],
['Axy','Dxz'],['Dxz','Byz'],['Byz','Axy'],
['Bxy','Bxz'],['Bxz','Ayz'],['Ayz','Bxy'],
['Bxy','Byz'],['Byz','Cxz'],['Cxz','Bxy'],
['Dxy','Axz'],['Axz','Dyz'],['Dyz','Dxy'],
['Dxy','Dxz'],['Dxz','Cyz'],['Cyz','Dxy'],
['Cxy','Bxz'],['Bxz','Dyz'],['Dyz','Cxy'],
['Cxy','Cyz'],['Cyz','Cxz'],['Cxz','Cxy']]

def mkrect(self,v1,v2,tag):
# add four rectangle defined by vectors to vertdict with
# keys 'A','B','C','D' + passed tag identifying plane
# of the rectangle i.e. 'xy','xz' or 'yz'

# Rectangle is: square when v1,v2 come from cuboctahedron
#    "       ": golden (1xphi) when v1,v2 from icosa
self.vertdict['A'+tag] =  v1 + v2
self.vertdict['B'+tag] = -v1 + v2
self.vertdict['C'+tag] = -v1 - v2
self.vertdict['D'+tag] =  v1 - v2

def drawedges(self,myfile):
# draw the edges defined in the edgepairs list
for i in self.edgepairs:
myfile.edge(self.vertdict[i[0]],self.vertdict[i[1]])

def drawrect(self, myfile, control=0):
# output the three rectangles (12 edges) to a povray file
# -- or filter to only draw edges with length self.control if
# control = 1 (used by Icosahedron drawself)
for tag in ["xy","xz","yz"]:
for labels in ['AB','BC','CD','DA']:
v1,v2 = self.vertdict[labels[0]+tag],self.vertdict[labels[1]+tag]
if control:
if length(v1-v2)==self.control:
myfile.edge(v1,v2)
else:
myfile.edge(v1,v2)

def volume(self):
# define a tetrahedron using one face and center
# and return 20 x its volume
v1 = Vector((0,0,0))
v2 = self.vertdict['Axy']
v3 = self.vertdict['Axz']
v4 = self.vertdict['Ayz']
oTetra = Tetra(v1,v2,v3,v4)
return 20.0 * oTetra.volume()

class Cubocta(Jitterbug,Shape):

def __init__(self):
# make 3 squares of root2 x root2 in the
# xy, xz and yz planes.  The resulting vertices
# will be those of the cuboctahedron

# xy plane
tag = "xy"
v1 = Vector((root2/2.0,      0.0, 0.0))
v2 = Vector((0.0      ,root2/2.0, 0.0))
self.mkrect(v1,v2,tag)

# xz plane
tag = "xz"
v1 = Vector((root2/2.0, 0.0, 0.0))
v2 = Vector((      0.0, 0.0, root2/2.0))
self.mkrect(v1,v2,tag)

# yz plane
tag = "yz"
v1 = Vector((0.0,       0.0, root2/2.0))
v2 = Vector((0.0, root2/2.0, 0.0))
self.mkrect(v1,v2,tag)

class Icosa(Jitterbug,Shape):

def __init__(self):
# make 3 golden rectangles of 1 x phi in the
# xy, xz and yz planes.  The resulting vertices
# will be those of the icosahedron

# xy plane
tag = "xy"
v1 = Vector((0.5,    0.0, 0.0))
v2 = Vector((0.0,phi/2.0, 0.0))
self.mkrect(v1,v2,tag)

# xz plane
tag = "xz"
v1 = Vector((phi/2.0, 0.0, 0.0))
v2 = Vector((    0.0, 0.0, 0.5))
self.mkrect(v1,v2,tag)

# yz plane
tag = "yz"
v1 = Vector((0.0, 0.0, phi/2.0))
v2 = Vector((0.0, 0.5, 0.0))
self.mkrect(v1,v2,tag)

def drawself(self,myfile):
self.drawedges(myfile)
self.drawrect(myfile,1)

class Octa(Jitterbug,Shape):

edgepairs = []  # don't need for this shape

def __init__(self):
# make 3 squares of 1 x 1 in the
# xy, xz and yz planes, with corners
# on the axis.  The resulting vertices
# will be those of the octahedron

# xy plane
tag = "xy"
v1 = Vector((root2/2.0,0.0,0.0))
v2 = Vector((0.0,root2/2.0,0.0))
self.mkrect(v1,v2,tag)

# xz plane
tag = "xz"
v1 = Vector((root2/2.0,0.0,0.0))
v2 = Vector((0.0,0.0,root2/2.0))
self.mkrect(v1,v2,tag)

# yz plane
tag = "yz"
v1 = Vector((0.0,root2/2.0,0.0))
v2 = Vector((0.0,0.0,root2/2.0))
self.mkrect(v1,v2,tag)

def mkrect(self,v1,v2,tag):
# add four rectangles defined by vectors to vertdict with
# keys 'A','B','C','D' + passed tag identifying plane
# of the rectangle i.e. 'xy','xz' or 'yz'

# Rectangle is: square when v1,v2 come from octahedron
# -- corners on the xyz axes

self.vertdict['A'+tag] =  v1
self.vertdict['B'+tag] =  v2
self.vertdict['C'+tag] = -v1
self.vertdict['D'+tag] = -v2

def drawself(self,myfile):
self.drawrect(myfile)

def volume(self):
# define a tetrahedron using one face and center
# and return 8 x its volume
v1 = Vector((0,0,0))
v2 = self.vertdict['Axy']
v3 = self.vertdict['Bxy']
v4 = self.vertdict['Byz']
oTetra = Tetra(v1,v2,v3,v4)
return 8.0 * oTetra.volume()

class Cube(Shape):

# Note: Duotet does not inherit from the Jitterbug class

def __init__(self):
self.vertdict = {'A': a,'B': b,'C': c,'D': d,
'E':-a,'F':-b,'G':-c,'H':-d}

self.edgepairs = [['A','F'],['A','G'],['A','H'],
['B','E'],['B','G'],['B','H'],
['C','E'],['C','F'],['C','H'],
['D','E'],['D','F'],['D','G']]

def volume(self):
v0,v1,v2,v3 =  self.getverts(['A','B','C','D'])
oTetra = Tetra(v0,v1,v2,v3)
return 3.0 * oTetra.volume()

class Tetra(Shape):

# Tetra may be defined using any four vectors.
# By default, it assumes the four vertices defined by
# closest packed unit-diameter spheres

# Note: Tetra does not inherit from the Jitterbug class

def __init__(self,v0=a,v1=b,v2=c,v3=d):
self.vertdict = {'A':v0,'B':v1,'C':v2,'D':v3}
self.edgepairs = [['A','B'],['A','C'],['A','D'],
['B','C'],['C','D'],['D','B']]

def volume(self):
# use scalar triple product (1/6)(a dot (b cross c)),
# with 8-folding because of unit = D (2 radii),
# and x 3rd power of synergetics constant to
# rationalize (thanks to Bucky Fuller)
v0,v1,v2,v3 =  self.getverts(['A','B','C','D'])
a = v0-v1
b = v0-v2
c = v0-v3
return (1.0/6.0)*abs(a.dot(b.cross(c)))*8*syn3

def volume1(self):
# evaluate above triple product using corresponding determinant
v0,v1,v2,v3 =  self.getverts(['A','B','C','D'])
e,x,y,z = [0]*3,[0]*3,[0]*3,[0]*3
e[0] = (v0-v1).xyz
e[1] = (v0-v2).xyz
e[2] = (v0-v3).xyz
for i in range(3):
x[i]=e[i][0]
y[i]=e[i][1]
z[i]=e[i][2]
det = (x[0]*y[1]*z[2]-x[0]*z[1]*y[2]-x[1]*y[0]*z[2]
+x[1]*z[0]*y[2]+x[2]*y[0]*z[1]-x[2]*z[0]*y[1])
return (1.0/6.0)*abs(det)*8*syn3

def volume2(self):

# a volume formula for the tetrahedron, using
# its six edge lengths as input (scalars, not vectors)

# thanks to Leonhard Euler and Gerald de Jong

# assuming sphere diameter as unit

v0,v1,v2,v3 =  self.getverts(['A','B','C','D'])
a = (v0-v1).length()
b = (v0-v2).length()
c = (v0-v3).length()
d = (v1-v2).length()
e = (v2-v3).length()
f = (v3-v1).length()

a2 = a*a
b2 = b*b
c2 = c*c
d2 = d*d
e2 = e*e
f2 = f*f

open   = ( f2 * a2 * b2
+ d2 * a2 * c2
+ a2 * b2 * e2
+ c2 * b2 * d2
+ e2 * c2 * a2
+ f2 * c2 * b2
+ e2 * d2 * a2
+ b2 * d2 * f2
+ b2 * e2 * f2
+ d2 * e2 * c2
+ a2 * f2 * e2
+ d2 * f2 * c2 )
closed = ( a2 * b2 * d2
+ d2 * e2 * f2
+ b2 * c2 * e2
+ a2 * c2 * f2 )
oppos  = ( a2 * e2 * (a2 + e2)
+ b2 * f2 * (b2 + f2)
+ c2 * d2 * (c2 + d2))

return (abs(open - closed - oppos)/2.0)**0.5

def volume3(self): # thanks to Tom Ace
v0,v1,v2,v3 =  self.getverts(['A','B','C','D'])
a0,a1,a2,a3 = (v0-v1).quadray()
b0,b1,b2,b3 = (v0-v2).quadray()
c0,c1,c2,c3 = (v0-v3).quadray()
det = ( a1*b2*c3 - a1*b3*c2 - b1*a2*c3 + b1*a3*c2
+ c1*a2*b3 - c1*a3*b2 - a0*b2*c3 + a0*b3*c2
+ a0*b1*c3 - a0*b1*c2 - a0*c1*b3 + a0*c1*b2
+ b0*a2*c3 - b0*a3*c2 - b0*a1*c3 + b0*a1*c2
+ b0*c1*a3 - b0*c1*a2 - c0*a2*b3 + c0*a3*b2
+ c0*a1*b3 - c0*a1*b2 - c0*b1*a3 + c0*b1*a2)
return (1.0/4.0)*abs(det)

# ------------------ print utilities -------

def qcoords(shape):
print " "
print "     quadray  (a,b,c,d)          "
print "---------------------------------"
for i in shape.vertdict.values():
print "(%f, %f, %f, %f)" % i.quadray()

def scoords(shape):
print " "
print "    spherical (r,phi,theta)      "
print "---------------------------------"
for i in shape.vertdict.values():
print "(%g, %g, %g)" % i.spherical()

def ccoords(shape):
print " "
print "    Cartesian (x,y,z)            "
print "---------------------------------"
for i in shape.vertdict.values():
print "(%f, %f, %f)" % (i.xyz)

```