```"""

A simple waterman maker for Jython

Finds all IVM balls within a radius such that volume is maximized, convexity preserved
(i.e. generates classic Waterman polyhedra).  Scales/colors for POV-Ray rendering.

Convex Hull finding added courtesy of John E. Lloyd, who implemented Qhull-style
algorithms in his QuickHull3D Java package -- and dependency for this program
(try same directory, as a folder named quickhull3d).
http://www.cs.ubc.ca/spider/lloyd/java/quickhull3d.html
Thanks also to Jim & Jython team, to Guido, Tim etc.

by Kirby Urner, June 16, 2005

Version 2.1

Note:  From my point of view, Steve Waterman is the originator of this family of
IVM-based polyhedra -- convex hulls with a maximum number of CCP ball vertices
of a given distance from the origin or less.  That's why I (without prompting
from Steve) named them Waterman polyhedra -- and the name appears to have stuck.

Usage:

Launch Jython
>>> import wpolys
>>> wpolys.waterman(10)
[ will write w10.pov to disk, ready for rendering with POV-Ray ]

I may not be able to provide any tech support.  For the hardy.

Note:  no set object in this version of Jython, so I use dictionaries
with None for values -- a way to preserve the uniqueness of keys.
"""
import quickhull3d as qh

def waterman(n):
"""
collect all ccp balls within a given sweepout radius
"""
filename = 'w'+str(n)+'.pov'
fileout  = open(filename,'w')

allballs = []
oldverts = []

# find balls from the outside in, 6 layers at a time,

# stop when convex hull registers no change in vertices

print 'Finding ivm lattice points'
layer = 0
for r in range(2*n,0,-2):
print '.',
balls = finder(r)
for b in balls:
newballs = permuter(*b)
allballs.extend( newballs )
layer += 1
if layer == 6 or r<=2:
hull = qh.QuickHull3D()
points = [qh.Point3d(*ball) for ball in allballs]
hull.build(points)
verts = [(p.x, p.y, p.z) for p in hull.getVertices()]
newverts = verts[:]
newverts.sort()
print '|',
if newverts == oldverts or r<=2:
break  # convex hull is stable and/or all points found

else:
oldverts = newverts
layer = 0

print "\nConvex hull found\n"
faces = [list(f) for f in hull.getFaces()]
edges = getedges(faces)

povarray(n, verts, fileout)

for face in faces:
povface(face, fileout)

for edge in edges:
povedge(edge, fileout)

fileout.close()

def finder(r):
"""
Explore first octant for unique balls
"""
# main loop

okballs = []
maxx = int(pow(r,0.5))
for x in range(maxx,0,-1):      # x drops from r to 1

for y in range(0,x+1):      # y up to x

if x**2 + y**2 > r:     # quit if too big

break
for z in range(0,y+1):  # z up to y

if not (x+y+z) % 2: # even sum only

if x**2 + y**2 + z**2 == r:
okballs.append((x,y,z)) # got one

return okballs

def permuter(x,y,z):
xyz = (x,y,z)
unique = {}  # import sets, use sets.Set() if < Python 2.4

# revert to dictionary for Jython (approx Python 2.1 as of now)

unique[(xyz[0],xyz[1],xyz[2])] = None
unique[(xyz[0],xyz[2],xyz[1])] = None
unique[(xyz[1],xyz[0],xyz[2])] = None
unique[(xyz[1],xyz[2],xyz[0])] = None
unique[(xyz[2],xyz[0],xyz[1])] = None
unique[(xyz[2],xyz[1],xyz[0])] = None
signs = [('+','+','-'),
('+','-','+'),
('-','+','+'),
('+','-','-'),
('-','+','-'),
('-','-','+'),
('-','-','-')]
for combo in unique.keys():
for pattern in signs:
newcombo = [eval(i + str(j))
for i,j in zip(pattern, combo)]
unique[tuple(newcombo)] = None
return unique.keys()

def getedges(faces):
alledges = {}
for f in faces:
for a,b in zip(f, f[1:]+[f[0]]):
candidate = [a,b]
candidate.sort()
alledges[tuple(candidate)] = None
return alledges.keys()

def povarray(n, verts, fileout, tex='texture6'):
"""
cut and paste into .pov file with standard (or custom) header
"""
scale = 1 # pow(2,0.5)*n

fileout.write(arrayskel % len(verts))
for b in verts[:-1]:
b = tuple([x/scale for x in b])
theline = "   <%d,%d,%d>,\n" % b
fileout.write(theline)
b = verts[-1]
b = tuple([x/scale for x in b])
fileout.write("   <%d,%d,%d>\n" % b)
fileout.write("}\n")
fileout.write( arrayloop % (0.04, 'texture1') )

def povface(face, fileout):
fileout.write("polygon{%s,\n" % (len(face)+1))
for p in face:
fileout.write("Pts[%s] " % p)
fileout.write("Pts[%s] " % face[0])
fileout.write("\ntexture {texture0}\n }\n")

def povedge(edge, fileout):
fileout.write("cylinder{ Pts[%s], Pts[%s], 0.04 pigment {texture1}}\n" % (edge[0], edge[1]))

//version 3.6
global_settings { assumed_gamma 2.2 }
#include "colors.inc"
#include "metals.inc"
#include "glass.inc"

#declare Cam_factor = %s;
#declare Camera_X = 0.5 * Cam_factor;
#declare Camera_Y = 0.5 * Cam_factor;
#declare Camera_Z = 1 * Cam_factor;

#declare texture0 = T_Copper_2A;
#declare texture1 = Col_Dark_Green;

camera { location  <Camera_X, Camera_Y, Camera_Z>
up        <0, 1.0,  0>    right     <-4/3, 0,  0>
direction <0, 0,  3>      look_at   <0, 0, 0>
rotate <0,0,0>}

light_source { <Camera_X - 2, Camera_Y + 5 , Camera_Z + 5> color White shadowless}
light_source { <Camera_X - 2, Camera_Y + 5 , Camera_Z - 3> color White shadowless}
light_source { <Camera_X + 2, Camera_Y - 5 , Camera_Z + 3> color White shadowless}

//Background:
background {color Black}

"""

arrayskel = """
# declare NumPoints = %s;
# declare Pts = array [NumPoints] {
"""

arrayloop = """
#declare ptno=0;
#while(ptno<NumPoints)
sphere{Pts[ptno] %s pigment{%s}}
#declare ptno = ptno+1;
#end
"""

```
`# code highlighted using py2html.py version 0.8`