By Kirby Urner, 4D Solutions, January, 2001

Note: coords.py and povray.py (used below) are in python101.zip
Other modules are part of the Standard Python Library.

import operator, os, string, povray, coords

# edit the line below to give proper command syntax on your computer
# comm0 = "g:\qhull\qhull < g:\python20\waterman.txt o  > qwater.txt"
# comm1 = "g:\qhull\qhull < g:\python20\waterman.txt FA > qwater.rep"

# or use qhull.bat in working subdirectory e.g. something like:
# g:\qhull\qhull < g:\python20\waterman.txt o  > qwater.txt
# g:\qhull\qhull < g:\python20\waterman.txt FA > qwater.rep
# pause
# and make comm0 = "qhull.bat"
# this is the approach taken below (comm1 commented out because not
# used) -- advantage is the inserted pause to let you read any DOS
# box messages (diagnostic step)

comm0 = "qhull.bat"

def mkwater(root,type=1):
    Generate a povray file for W(root) -- this is the main program
    If type is not 'w', then fixed radius option is invoked
    allcoords = []
    poscoords = poly(root)    

    if len(poscoords)==0:
        print "Missing for root %s\n" % root
        return 0

    if type==1:
        for i in range(1,10):
            if root>(i+1):
                poscoords += poly(root-i)        

    for v in poscoords:        
        allcoords = allcoords + negperm(v)
    print "Raw data generated for root %s" % root

    print "Convex hull generated"
    # os.system(comm1)
    print "Summary statistics filed"

    verts,faces = readfaces()
    nbverts,nbfaces,volume = readrep()

    print ("Vertices: %s  Faces: %s  Edges: %s  Volume: %s" %
          (nbverts,nbfaces,nbverts + nbfaces - 2, volume * 3))

    print "Povray file generated\n"        
    return 1

def listbase(root):
    scalefactor = 1/(2.0)**0.5    
    print "Distance = Root(%s)" % root
    allcoords = []
    poscoords = poly(root)
    for v in poscoords:        
        allcoords = allcoords + negperm(v)
    vectors = map(coords.Vector,allcoords)  # make vector objects
    for v in vectors:
        base = (v.rotz(45.0))*scalefactor            
        print "(%g,%g,%g)" % (base.xyz[0],base.xyz[1],base.xyz[2])
def poly(root,perm=1):
    List all-positive, integral (x,y,z) positions at the centers
    of CCP spheres that are distance pow(2*root,0.5) from the
    origin -- optionally suppress permutations with 2nd argument
    of 0.

    CCP spheres have diameter root(2) in this coordinate system,
    so 2nd powers of sphere-center distances from (0,0,0) will
    be multiples of 2, i.e. pow(N*root(2),2) = 2*(N**2).  N**2
    = root, the first argument to this method.  A scale factor
    of 1/root(2) may be applied later such that root = 3 results
    in CCP spheres a distance square root of 3 from (0,0,0).

    x + y + z = even number (i.e. (x+y+z)%2==0.  This stipulation
    ensures that our integer (x,y,z) coordinates are coincident
    with the center of a CCP sphere.

    See: http://www.inetarena.com/~pdx4d/ocn/wmalgorithm.html for
    additional background information
    rad = root*2
    coords = []
    for x in range(int(pow(rad,.5)),-1,-1):  # descending order
        if x**2 < rad/3:
        for y in range(x+1):
            if x**2 + y**2 > rad:
            for z in range(y+1):
                if (x+y+z)%2 == 0:  # check for even sum
                    if x**2 + y**2 + z**2 > rad:
                    if x**2 + y**2 + z**2 == rad:
    if len(coords)==0: return coords                    
    elif perm: return reduce(operator.add,map(permute,coords))
    else:      return coords

def permute(e):
    return list of unique permutations of 3 elements
    result = []    
    for i in range(3):
        for j in range(3):
            for k in range(3):
                if i<>j and j<>k and i<>k:
                   permutation = (e[i],e[j],e[k])
                   if not permutation in result:
    return result

def negperm(i):
    return list of tuples permuted with negative signs
    def addperm(r,i):
        if not i in r:
        return r
    r = [i]
    r = addperm(r,( i[0],-i[1], i[2]))
    r = addperm(r,( i[0], i[1],-i[2]))
    r = addperm(r,( i[0],-i[1],-i[2]))    
    r = addperm(r,(-i[0], i[1], i[2]))
    r = addperm(r,(-i[0],-i[1], i[2]))
    r = addperm(r,(-i[0], i[1],-i[2]))
    r = addperm(r,(-i[0],-i[1],-i[2]))

    return r

def writepoints(verts):
    Write the vertices, computed and permuted, to a text
    file in a format suitable for processing by Qhull
    txtfileobj = open("waterman.txt",'w')

    txtfileobj.write("3 #3D file\n")

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

def readfaces():
    Read the data file output by Qhull, analyzing it into
    a list of vertices and a list of faces
    txtfileobj = open("qwater.txt",'r') # qhull output file    
    datalines = txtfileobj.readlines()
    nbverts,nbfaces,nbedges = map(int,(string.split(datalines[1])))
    verts,faces = [],[]
    for i in range(2,nbverts+2):

    for i in range(nbverts+2,nbverts+nbfaces+2):

    return verts, faces

def readrep():
    Read the data file output by Qhull, analyzing it into
    a list of vertices and a list of faces
    txtfileobj = open("qwater.rep",'r') # qhull output file    
    datalines = txtfileobj.readlines()
    nbverts = float((string.split(datalines[3])[3]))
    nbfaces = float((string.split(datalines[4])[3]))
    volume  = float((string.split(datalines[15])[2]))


    return nbverts,nbfaces,volume

def genpovray(verts, faces, root, type):
    Generate the actual Povray file, scaled by root (note:
    all computations have been with integers up to this point)
    scalefactor = 1/(2.0*root)**0.5 
    vectors = map(coords.Vector,verts)  # make vector objects
    if type==1:

        print "Color coding vertices"
        for i in range(len(vectors)):
           if vectors[i].dot(vectors[i]) == 2*root:
               vectors[i].color = "Gold"
           elif vectors[i].dot(vectors[i]) == 2*(root-1):
               vectors[i].color = "Silver"
           elif vectors[i].dot(vectors[i]) == 2*(root-2):
               vectors[i].color = "Green"
           elif vectors[i].dot(vectors[i]) == 2*(root-3):
               vectors[i].color = "Blue"
           elif vectors[i].dot(vectors[i]) == 2*(root-4):
               vectors[i].color = "Yellow"
           elif vectors[i].dot(vectors[i]) == 2*(root-5):
               vectors[i].color = "Magenta"
           elif vectors[i].dot(vectors[i]) == 2*(root-6):
               vectors[i].color = "LightSteelBlue"
           elif vectors[i].dot(vectors[i]) == 2*(root-7):
               vectors[i].color = "Pink"
           elif vectors[i].dot(vectors[i]) == 2*(root-8):
               vectors[i].color = "DarkTan"
           elif vectors[i].dot(vectors[i]) == 2*(root-9):
               vectors[i].color = "DustyRose"
               vectors[i].color = "Black"
        myfile = povray.Povray("wm"+str(root)+".pov",bgcolor="Black")
        myfile = povray.Povray("fr"+str(root)+".pov")
    newvectors = map(operator.mul,len(vectors)*[scalefactor],vectors)

    if type==1:
        for i in range(len(newvectors)):        
           newvectors[i].drawn = 0
           newvectors[i].color = vectors[i].color
    vectors = newvectors    
    myfile.facecolor = "Cyan"
    myfile.cylcolor  = "Red"
    myfile.sphradius = 0.05
    myfile.cylradius = 0.02
    for face in faces:
        vlist = []
        prev = 0
        for v in face[1:]:
            if prev<>0:
            prev = vectors[v]
            if type==1 and prev.drawn==0:
               myfile.sphcolor = prev.color
               prev.drawn = 1