```"""
VPYTHON VERSION

Originally for POV-Ray and VRML work.  Repurposed to work
directly with VPython.

By Kirby Urner, Oregon Curriculum Network

Ver 1.4:  March 15, 2005

Ver 0.2:  May 29, 2000

"""

from visual import *
from operator import mul, sub
from math import sqrt,pi,sin,cos
import os, string, copy

#### Python 2.3 fixes -- added by JMZ ####

try:
set
except:
from sets import Set as set

try:
sorted
except:
def sorted(lst):
lst2 = list(lst)
lst2.sort()
return lst2

##########################################

colordict = {"Red"    :"1.0 0.0 0.0",
"Green"  :"0.0 1.0 0.0",
"Blue"   :"0.0 0.0 1.0",
"Yellow" :"1.0 1.0 0.0",
"Cyan"   :"0.0 1.0 1.0",
"Magenta":"1.0 0.0 1.0",
"White"  :"1.0 1.0 1.0",
"Black"  :"0.0 0.0 0.0",
"Orange" :"1.0 0.5 0.0",
"Violet" :"0.309804 0.184314 0.309804",
"Indigo" :"0.294117 0.0      0.509803",
"Brown"  :"0.647059 0.164706 0.164706",
"Wood"   :"0.647059 0.164706 0.164706",
"Grey"   :"0.752941 0.752941 0.752941",
"Gray"   :"0.752941 0.752941 0.752941",
"Gold"   :"0.8 0.498039 0.196078",
"Silver" :"0.90 0.91 0.98"}

root2   = 2.0**0.5

class Qvector(object):
"""
Converts Quadray coordinates to XYZ (from coords.py)

Used to be a subclass of a fully functional vector but
here is just used to convert from one coordinate system
to another
"""

def __init__(self,arg):

self.coords = self.norm(arg)

a,b,c,d     =  self.coords
self.xyz    = ((0.5/root2) * (a - b - c + d),
(0.5/root2) * (a - b + c - d),
(0.5/root2) * (a + b - c - d))

def norm(self,plist):
"""Normalize such that 4-tuple all non-negative members."""
return tuple(map(sub,plist,[min(plist)]*4))

#==================[ Points of Interest ]====================

"""
* 26 data points A-Z define six polyhedra in the concentric hierarchy

Labels of   Numbers of
Shape              Volume  Vertices    Vertices, Edges, Faces
-------------------------------------------------------------
Tetrahedron           1      A -D          4       6      4
Inv Tetra             1      E -H          4       6      4
Duo-tet Cube          3      A -H          8      12      6
Octahedron            4      I -N          6      12      8
Rh Dodecahedron       6      A -N         14      24     12
Cuboctahedron        20      O -Z         12      24     14

See: http://www.inetarena.com/~pdx4d/ocn/graphics/povlabels.gif

Using the quadray apparatus (4 basis vectors to the corners
of a regular tetrahedron with edges = 1 sphere diameter)

"""

Vset = {}

def init():
global Vset

ORIGIN = vector(0,0,0)

A = vector(Qvector((1,0,0,0)).xyz)  # center to corner of tetrahedron

B = vector(Qvector((0,1,0,0)).xyz)  #          "

C = vector(Qvector((0,0,1,0)).xyz)  #          "

D = vector(Qvector((0,0,0,1)).xyz)  #          "

# tetrahedron's dual (also a tetrahedron i.e. inverted tet)

E,F,G,H = B+C+D,A+C+D,A+B+D,A+B+C

# tetrahedron + dual (inverted tet) = duo-tet cube

# octahedron vertices from pairs of tetrahedron radials

I,J,K,L,M,N = A+B, A+C, A+D, B+C, B+D, C+D

# octahedron + dual (cube) = rhombic dodecahedron

# cuboctahedron vertices from pairs of octahedron radials

O,P,Q,R,S,T = I+J, I+K, I+L, I+M, N+J, N+K
U,V,W,X,Y,Z = N+L, N+M, J+L, L+M, M+K, K+J

"""
* the Jitterbug Transformation creates a basis for additional vertices

Labels of    Numbers of
Shape              Volume  Vertices     Vertices, Edges, Faces
---------------------------------------------------------------
Rh Triacontahedron    5*     a -t,O1-Z1   32       60     30
Reg Dodecahedron    ~15.21   a -t         20       30     12
Icosahedron         ~18.51   O1-Z1        12       30     20

* initialization includes scaling/shrinking
"""

# icosahedron vertices from jitterbugged cubocta vertices

phi   = (1.+5**0.5)/2.
k1    = phi/2.
k2    = 0.5
k3    = (3.**0.5)/2.
ef    = 1.0/phi
tau   = ef
tf    = pow( (sqrt(2) * (2.0+sqrt(5))/6.0),(1.0/3.0))

def jitterbug((v1,v2)): return k1*norm(v1) + k2*norm(v2)

O1,P1,Q1,R1,S1,T1 = map(jitterbug,[(I,J),(K,I),(L,I),(I,M),(N,J),(K,N)])
U1,V1,W1,X1,Y1,Z1 = map(jitterbug,[(L,N),(N,M),(J,L),(M,L),(M,K),(J,K)])

# pentagonal dodeca vertices from icosa triplets

def dualicosa((v1,v2,v3)): return k3*norm((v1+v2+v3))

a,b,c,d = map(dualicosa,([O1,Q1,W1],[O1,Z1,P1],[P1,R1,Y1],[R1,Q1,X1]))
e,f,g,h = map(dualicosa,([W1,S1,U1],[U1,X1,V1],[Y1,T1,V1],[Z1,T1,S1]))
i,j,k,l = map(dualicosa,([O1,P1,R1],[O1,Q1,R1],[O1,Z1,W1],[W1,S1,Z1]))
m,n,o,p = map(dualicosa,([Q1,W1,U1],[Q1,U1,X1],[Z1,P1,T1],[Y1,P1,T1]))
q,r,s,t = map(dualicosa,([R1,Y1,X1],[V1,Y1,X1],[T1,S1,V1],[U1,S1,V1]))

verts = locals()
for v in verts:
Vset[v] = verts[v]

init()

#====================[ Polyhedra ]====================

class Shape:

sphcolor  = "Red"
cylcolor  = "Orange"
facecolor = "Cyan"

cyltexture = ""
facetexture = ""
sphtexture = ""

showfaces = 0
showedges = 1
showverts = 1
showlabels = 0

def __init__(self):
self.getedges()
self.getvertices()

def getvertices(self):
"""Extract vertices from edges list.

Add unique string:vector pairs to the vertices
dictionary by stepping through edge list, and
stop any time the total number of vertices is reached
"""
self.vertices = {} # initialize dictionary

nvertices =  len(self.edges) - len(self.faces) + 2 # V = E - F + 2

for edge in self.edges:      # e.g. edge = ['A','B']

for label in edge:       # e.g. j = 'A'

v = Vset[label] # initialize vectors using V.label

self.vertices[label] = vector(v.x, v.y, v.z) # eval(label)  # vector <-- eval(string)

if len(self.vertices)==nvertices:
break
self.center = vector(0,0,0)

def getedges(self):
"""
Extract edges from the faces list.
"""
edges = set()
for f in self.faces:
pairs = zip(f , f[1:]+(f[0],))
for p in pairs:
self.edges = sorted(list(edges))

def getangles(self):
"""
List surface angles as vertex 3-tuples
"""
anglelist = []
for face in self.faces:
for j in range(len(face)):
angle = (face[j],face[j-1],face[j-2])
anglelist.append(angle)
return anglelist

def angles(self):
"""
List surface angles in degrees
"""
angles = []
for angle in self.getangles():
v1 = self.vertices[angle[2]] - self.vertices[angle[1]]
v2 = self.vertices[angle[0]] - self.vertices[angle[1]]
angles.append(v1.diff_angle(v2))
return angles

def scale(self,factor):
"""
Resize shape by factor.

Translate to the origin, scale, translate back.
Multiply volume by 3rd power of scale factor.
"""
for label in self.vertices.keys():
self.vertices[label] = self.vertices[label] * factor
self.volume = self.volume * pow(factor,3)
return self

__mul__  = scale
__rmul__ = __mul__

def rotate(self,axis,degrees):
"""
Rotate shape in place around axis X Y or Z by some angle

Rotation matrices specified below.  Translate to the
origin, rotate, translate back
"""
newself = copy.deepcopy(self)
current = newself.center
newself.move(-current)  # go to origin

for label in newself.vertices.keys():
if   axis == 'X':
newself.vertices[label]= rotx(newself.vertices[label],theta)
elif axis == 'Y':
newself.vertices[label]= roty(newself.vertices[label],theta)
elif axis == 'Z':
newself.vertices[label]= rotz(newself.vertices[label],theta)

return newself

def move(self,v):
"""
Add vector v to each vertex
"""
for label in self.vertices.keys():
self.vertices[label] = self.vertices[label] + v
self.center = self.center + v

def translate(self,v):
"""
Translate shape by vector v.
Return new shape.
"""
newself = copy.deepcopy(self)
newself.move(v)
return newself

def draw(self, trace=False):
"""
Output polyhedron to VPython
"""
self.garbage = []
self.cyls = {}
self.verts = {}

ecolor = tuple([float(i) for i in colordict[self.cylcolor].split()])
vcolor = tuple([float(i) for i in colordict[self.sphcolor].split()])

for a,b in self.edges:
v0 = self.vertices[a]
v1 = self.vertices[b]
self.garbage.append(sph)
self.verts[(a,b)] = sph
if trace:
else:
self.garbage.append(cyl)
self.cyls[(a,b)]=cyl
self.verts[(b,a)] = sph
self.garbage.append(sph)

def redraw(self):
for a,b in self.cyls:
self.cyls[(a,b)].pos  = self.vertices[a]
self.cyls[(a,b)].axis = self.vertices[b] - self.vertices[a]
self.verts[(a,b)].pos = self.vertices[a]
self.verts[(b,a)].pos = self.vertices[b]

def _trace(self, a, b, r, c, therate=500):
v0 = self.vertices[a]
v1 = self.vertices[b]
start = cylinder(pos=v0, length=0, color=c, radius = r)
self.garbage.append(start)
self.cyls[(a,b)] = start
v3 = v1-v0
for i in range(1,101):
rate(therate)
start.axis = v3*i/100.0

def delete(self,therate=200):
for obj in self.garbage:
if not therate==0:
rate(therate)
obj.visible = 0

class XYZaxes(Shape):
edges = [('ORIGIN','I'),('ORIGIN','J'),('ORIGIN','K'),
('ORIGIN','L'),('ORIGIN','M'),('ORIGIN','N')]

sphcolor = cylcolor = "Green"

faces = []

def __init__(self):
self.getvertices()
self.volume  = 0.0

class Qaxes(Shape):

edges = [('ORIGIN','A'),('ORIGIN','B'),
('ORIGIN','C'),('ORIGIN','D')]

sphcolor = cylcolor = "Yellow"

faces = []

def __init__(self):
self.getvertices()
self.volume   = 0.0

class Invqaxes(Shape):

edges = [('ORIGIN','E'),('ORIGIN','F'),
('ORIGIN','G'),('ORIGIN','H')]

faces = []

def __init__(self):
self.getvertices()
self.volume   = 0.0

class Tetra(Shape):
"""
Labels of   Numbers of
Shape          Volume  Vertices    Vertices, Edges, Faces
---------------------------------------------------------
Tetrahedron        1      A-D           4       6      4
"""

faces = [('A','B','C'),('A','C','D'),('A','D','B'),('B','C','D')]

sphcolor = cylcolor = "Orange"

def __init__(self):
Shape.__init__(self)
self.volume   = 1.0

class Invtetra(Shape):
"""
Labels of   Numbers of
Shape          Volume  Vertices    Vertices, Edges, Faces
---------------------------------------------------------
Inv Tetra          1      E-H           4       6      4
"""

faces = [('E','F','G'),('E','G','H'),('E','H','G'),('F','G','H')]

sphcolor = cylcolor = "Black"

def __init__(self):
Shape.__init__(self)
self.volume   = 1.0

class Cube(Shape):
"""
Labels of   Numbers of
Shape          Volume  Vertices    Vertices, Edges, Faces
---------------------------------------------------------
Duo-tet Cube       3      A-H           8      12      6
"""

faces = [('A','H','C','F'),('A','H','B','G'),('B','E','C','H'),
('B','E','D','G'),('D','G','A','F'),('C','E','D','F')]

sphcolor = cylcolor = "Green"

def __init__(self):
Shape.__init__(self)
self.volume   = 3.0

class DST_Cube1(Cube):

"""
Used with Design Science Toys pix
"""

showedges = 1
showverts = 1
showfaces = 0
sphtexture = "T_Chrome_2C"
cyltexture = "T_Grnt15"

class Octa(Shape):
"""
Labels of   Numbers of
Shape          Volume  Vertices    Vertices, Edges, Faces
---------------------------------------------------------
Octahedron         4      I-N           6      12      8
"""

faces = [('I','J','K'),('I','L','M'),('I','J','L'),('I','M','K'),
('N','J','K'),('N','K','M'),('N','M','L'),('N','L','J')]

sphcolor = cylcolor = "Red"

def __init__(self):
Shape.__init__(self)
self.volume   = 4.0

class Rhdodeca(Shape):
"""
Labels of   Numbers of
Shape          Volume  Vertices    Vertices, Edges, Faces
---------------------------------------------------------
Rh Dodecahedron    6      A-N          14      24     12
"""

faces = [('N','F','J','C'),('N','F','K','D'),('N','D','M','E'),
('N','C','L','E'),('I','A','J','H'),('I','H','L','B'),
('I','B','M','G'),('I','G','K','A'),('K','F','J','A'),
('K','D','M','G'),('M','B','L','E'),('L','C','J','H')]

sphcolor = cylcolor = "Blue"

def __init__(self):
Shape.__init__(self)
self.volume   = 6.0

class DST_Rhdodeca(Rhdodeca):
"""
Used with Design Science Toys pix
"""

showedges = 0
showverts = 0
showfaces = 1
facetexture = "T_Wood25"

class DST_Rhdodeca1(Rhdodeca):

"""
Used with Design Science Toys pix
"""

showedges = 1
showverts = 1
showfaces = 0
sphtexture = "T_Chrome_2C"
cyltexture = "T_Grnt15"

class Cubocta(Shape):
"""
Labels of   Numbers of
Shape          Volume  Vertices    Vertices, Edges, Faces
---------------------------------------------------------
Cuboctahedron     20      O-Z          12      24     14
"""

sphcolor = cylcolor = "Yellow"

faces = [('O','Q','W'),('O','Z','P'),('P','R','Y'),('R','Q','X'),
('W','S','U'),('U','X','V'),('Y','T','V'),('Z','T','S'),
('O','P','R','Q'),('O','W','S','Z'),('Q','W','U','X'),
('P','Z','T','Y'),('R','X','V','Y'),('T','S','U','V')]

def __init__(self):
Shape.__init__(self)
self.volume   = 20.0

class VEspokes(Shape):

sphcolor = cylcolor = "Violet"

edges = [('ORIGIN','O'),('ORIGIN','P'),
('ORIGIN','Q'),('ORIGIN','R'),
('ORIGIN','S'),('ORIGIN','T'),
('ORIGIN','U'),('ORIGIN','V'),
('ORIGIN','W'),('ORIGIN','X'),
('ORIGIN','Y'),('ORIGIN','Z')]

faces = []

def __init__(self):
self.getvertices()
self.volume   = 0.0

class Icosa(Shape):
"""
Labels of   Numbers of
Shape          Volume  Vertices    Vertices, Edges, Faces
---------------------------------------------------------
Icosahedron    ~18.51    O1-Z1        12      30     20
"""

sphcolor = cylcolor = "Cyan"

faces = [('O1','Q1','W1'),('O1','Z1','P1'),('P1','R1','Y1'),('R1','Q1','X1'),
('W1','S1','U1'),('U1','X1','V1'),('Y1','T1','V1'),('Z1','T1','S1'),
('O1','P1','R1'),('O1','Q1','R1'),('O1','Z1','W1'),('W1','S1','Z1'),
('Q1','W1','U1'),('Q1','U1','X1'),('Z1','P1','T1'),('Y1','P1','T1'),
('R1','Y1','X1'),('V1','Y1','X1'),('T1','S1','V1'),('U1','S1','V1')]

def __init__(self):
Shape.__init__(self)
self.volume   = 0

class Jiticosa(Shape):
"""
Labels of   Numbers of
Shape          Volume  Vertices    Vertices, Edges, Faces
---------------------------------------------------------
Icosahedron    ~18.51    O1-Z1        12      30     20
"""

sphcolor = cylcolor = "Cyan"

faces = [('O1','Q1','W1'),('O1','Z1','P1'),('P1','R1','Y1'),('R1','Q1','X1'),
('W1','S1','U1'),('U1','X1','V1'),('Y1','T1','V1'),('Z1','T1','S1'),
('O1','P1','R1','Q1'),('O1','W1','S1','Z1'),('Q1','W1','U1','X1'),
('P1','Z1','T1','Y1'),('R1','X1','V1','Y1'),('T1','S1','U1','V1')]

def __init__(self):
Shape.__init__(self)
self.volume   = 0

class Icospokes(Shape):

sphcolor = cylcolor = "Yellow"

edges = [('ORIGIN','O1'),('ORIGIN','P1'),
('ORIGIN','Q1'),('ORIGIN','R1'),
('ORIGIN','S1'),('ORIGIN','T1'),
('ORIGIN','U1'),('ORIGIN','V1'),
('ORIGIN','W1'),('ORIGIN','X1'),
('ORIGIN','Y1'),('ORIGIN','Z1')]

faces = []

def __init__(self):
self.getvertices()
self.volume   = 0.0

class Grects(Shape):

faces = [('Q1','P1','T1','U1'),('X1','Y1','Z1','W1'),('S1','V1','R1','O1')]
sphcolor = cylcolor = "Violet"

def __init__(self):
Shape.__init__(self)
keylen = self.length(('P1','Q1'))
self.edges = [e for e in self.edges if not self.length(e)<keylen ]
self.volume   = 0.0

def length(self,e):
return length(self.vertices[e[0]]-self.vertices[e[1]])

class Pentdodeca(Shape):

faces = [('o','p','c','i','b'),('r','f','n','d','q'),('n','f','t','e','m'),
('e','l','k','a','m'),('j','a','k','b','i'),('q','d','j','i','c'),
('j','d','n','m','a'),('p','g','r','q','c'),('g','s','t','f','r'),
('k','l','h','o','b'),('h','l','e','t','s'),('h','s','g','p','o')]

sphcolor = cylcolor = "Brown"

def __init__(self):
Shape.__init__(self)
self.volume   = 0

class Pentcube(Shape):

sphcolor = cylcolor = "Green"

faces = [('j','b','p','q'),('j','q','f','m',),('j','b','l','m'),('l','s','p','b'),
('f','s','p','q'),('f','s','l','m')]

def __init__(self):
Shape.__init__(self)
self.volume = 0

class Penttet(Shape):

sphcolor = cylcolor = "Green"

faces = [('q','b','s'),('s','q','m'),('m','b','q'),('m','b','s')]

def __init__(self):
Shape.__init__(self)
self.volume = 0

class Rhtriaconta(Shape):

faces = [('o', 'P1', 'b', 'Z1'),('i', 'O1', 'b', 'P1'),
('O1', 'k', 'Z1', 'b'),('e', 'U1', 't', 'S1'),
('W1', 'm', 'U1', 'e'),('l', 'W1', 'e', 'S1'),
('Z1', 'l', 'S1', 'h'),('W1', 'l', 'Z1', 'k'),
('c', 'P1', 'p', 'Y1'),('c', 'R1', 'i', 'P1'),
('Y1', 'q', 'R1', 'c'),('a', 'W1', 'k', 'O1'),
('q', 'X1', 'd', 'R1'),('r', 'X1', 'q', 'Y1'),
('n', 'X1', 'f', 'U1'),('R1', 'j', 'O1', 'i'),
('Q1', 'd', 'X1', 'n'),('Q1', 'n', 'U1', 'm'),
('a', 'Q1', 'm', 'W1'),('Q1', 'j', 'R1', 'd'),
('j', 'Q1', 'a', 'O1'),('r', 'V1', 'f', 'X1'),
('f', 'V1', 't', 'U1'),('V1', 's', 'S1', 't'),
('g', 'V1', 'r', 'Y1'),('T1', 'p', 'P1', 'o'),
('h', 'T1', 'o', 'Z1'),('s', 'T1', 'h', 'S1'),
('T1', 's', 'V1', 'g'),('T1', 'g', 'Y1', 'p')]

sphcolor = cylcolor = "Magenta"

def __init__(self):
Shape.__init__(self)
self.volume   = 0

class Pentagon(Shape):

global p1,p2,p3,p4,p5

p1 = vector([cos54,  sin54, 0])
p2 = vector([cos126, sin126,0])
p3 = vector([cos198, sin198,0])
p4 = vector([0,-1,0])
p5 = vector([cos342, sin342,0])

faces = [('p1','p2','p3','p4','p5'),('p1','p3'),('p1','p4')]

def __init__(self):
Shape.__init__(self)
self.volume   = 0

class Pentagon2(Pentagon):

faces = [('p1','p2','p3','p4','p5')]

def __init__(self):
Pentagon.__init__(self)
self.volume   = 0

class Square(Shape):

global s1,s2,s3,s4

s1 = vector([ sin45,  sin45, 0])
s2 = vector([ sin45, -cos45, 0])
s3 = vector([-sin45, -cos45, 0])
s4 = vector([-sin45,  sin45, 0])

faces = [('s1','s2','s3','s4'),('s1','s3')]

def __init__(self):
Shape.__init__(self)
self.volume = 0

class Square2(Square):

faces = [('s1','s2','s3','s4')]

def __init__(self):
Square.__init__(self)
self.volume = 0

#=====================[ Matrix Ops ]==========================

def rotx(v,theta):
"""
Multiply v by X-axis rotation matrix
"""
return vector(v, theta, axis=(1,0,0))

def roty(v,theta):
"""
Multiply v by Y-axis rotation matrix
"""
return vector(v, theta, axis=(0,1,0))

def rotz(v,theta):
"""
Multiply v by Z-axis rotation matrix
"""
return vector(v, theta, axis=(0,0,1))

def det(qlist):
a0,a1,a2,a3 = qlist[0].coords
b0,b1,b2,b3 = qlist[1].coords
c0,c1,c2,c3 = qlist[2].coords
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 det

def util():
"""
Write the vertices, computed and permuted, to a text
file in a format suitable for processing by Qhull
"""

verts = [a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,
O1,P1,Q1,R1,S1,T1,U1,V1,W1,X1,Y1,Z1]

txtfileobj = open("tria.txt",'w')

txtfileobj.write("3 #3D file\n")
txtfileobj.write(str(len(verts))+"\n")

for vert in verts:
txtfileobj.write("%r %r %r\n" % vert.xyz)

txtfileobj.close()

comm0 = '"d:\program files\qhull\qhull" < "d:\program files\python21\tria.txt" o  > qtria.txt'
os.system(comm0)

"""
Read the data file output by Qhull, analyzing it into
a list of vertices and a list of faces
"""

verts = []
faces = []

txtfileobj = open(filename,'r') # qhull output file

`# code highlighted using py2html.py version 0.8`