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

writepoints(allcoords)
print "Raw data generated for root %s" % root

os.system(comm0)
print "Convex hull generated"

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

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

genpovray(verts,faces,root,type)
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
"""
coords = []
for x in range(int(pow(rad,.5)),-1,-1):  # descending order
break
for y in range(x+1):
if x**2 + y**2 > rad:
break
for z in range(y+1):
if (x+y+z)%2 == 0:  # check for even sum
if x**2 + y**2 + z**2 > rad:
break
if x**2 + y**2 + z**2 == rad:
coords.append((x,y,z))

if len(coords)==0: return 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:
result.append(permutation)
return result

def negperm(i):
"""
return list of tuples permuted with negative signs
"""
if not i in r:
r.append(i)
return r

r = [i]

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")
txtfileobj.write(str(len(verts))+"\n")

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

txtfileobj.close()

"""
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
nbverts,nbfaces,nbedges = map(int,(string.split(datalines[1])))

verts,faces = [],[]

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

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

txtfileobj.close()
return verts, faces

"""
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
nbverts = float((string.split(datalines[3])[3]))
nbfaces = float((string.split(datalines[4])[3]))
volume  = float((string.split(datalines[15])[2]))

txtfileobj.close()

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"
else:
vectors[i].color = "Black"

myfile = povray.Povray("wm"+str(root)+".pov",bgcolor="Black")

else:
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"

for face in faces:
vlist = []
prev = 0
for v in face[1:]:
vlist.append(vectors[v])
if prev<>0:
myfile.edge(prev,vectors[v])
prev = vectors[v]
if type==1 and prev.drawn==0:
myfile.sphcolor = prev.color
prev.drawn = 1
myfile.point(prev)
myfile.edge(vlist[-1],vlist[0])
myfile.face(vlist)

myfile.close()
```