"""

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.  
Use at your own risk.

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')
    fileout.write(header % (10 * pow(n,0.5)))
    
    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]))
    
header= """//POV-Ray script
//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
"""