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