"""
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
Last modified:  March 14, 2005
Last modified:  April 15, 2001

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"}

deg2rad = math.pi/180
rad2deg = 180/math.pi
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)

See:  http://www.teleport.com/~pdx4d/quadray/quadray.gif
"""

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:

    sphradius = 0.02
    cylradius = 0.02

    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:
	        edges.add(tuple(sorted(p)))
	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
        theta = degrees*deg2rad           
        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)
            newself.move(current)   # return to current position
        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]            
            sph = sphere(pos=v0, radius=self.sphradius, color = vcolor)            
            self.garbage.append(sph)
            self.verts[(a,b)] = sph
            if trace:
                self._trace(a, b, self.cylradius, ecolor)
            else:
                cyl=cylinder(pos=v0, axis=v1-v0, radius=self.cylradius, color = ecolor)
                self.garbage.append(cyl)
                self.cyls[(a,b)]=cyl
            sph = sphere(pos=v1, radius=self.sphradius, color = vcolor)
            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"
    sphradius = 0.05
    
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"
    sphradius = 0.05
    
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   = 20*vol1((O1,Q1,W1))

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   = 20*vol1((O1,Q1,W1))
        
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   = 24. * vol1((h,l,e))  + 12. * vol1((h,e,t))

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   = 60. * vol1((o,P1,Z1))
        
class Pentagon(Shape):

    global p1,p2,p3,p4,p5

    sin54,  cos54  = sin( 54*deg2rad), cos( 54*deg2rad)
    sin126, cos126 = sin(126*deg2rad), cos(126*deg2rad)
    sin198, cos198 = sin(198*deg2rad), cos(198*deg2rad)
    sin342, cos342 = sin(342*deg2rad), cos(342*deg2rad)        

    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

    sin45,cos45 = sin(45*deg2rad),cos(45*deg2rad)

    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 ]==========================
        
deg2rad = pi/180

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)
    
def readfaces(filename="qtria.txt"):
    """
    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    
    datalines = txtfileobj.readlines()
    nbverts,nbfaces,nbedges = map(int,(string.split(datalines[1])))
        
    for i in range(2,nbverts+2):
        verts.append(map(float,string.split(datalines[i])))

    for i in range(nbverts+2,nbverts+nbfaces+2):
        faces.append(map(int,string.split(datalines[i]))[1:])

    txtfileobj.close()
    return verts, faces