"""
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"
verts,faces = readfaces()
nbverts,nbfaces,volume = readrep()
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
additional background information
"""
rad = root*2
coords = []
for x in range(int(pow(rad,.5)),-1,-1): # descending order
if x**2 < rad/3:
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
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:
result.append(permutation)
return result
def negperm(i):
"""
return list of tuples permuted with negative signs
"""
def addperm(r,i):
if not i in r:
r.append(i)
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")
txtfileobj.write(str(len(verts))+"\n")
for vert in verts:
txtfileobj.write("%s %s %s\n" % vert)
txtfileobj.close()
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):
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
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]))
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"
myfile.sphradius = 0.05
myfile.cylradius = 0.02
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()