Hi Steve,
Thanks. From
>>>mkpoints.py 32 32.txt
here's a sample of the output
-0.396087591388 -0.781284022758 0.482400140683
-0.967387012461 -0.0838047084421 -0.239037944614
0.0208969821213 -0.489420208746 0.871797668848
0.887250003871 -0.258893773768 -0.38178717178
0.426352071227 -0.457758408728 -0.780180203927
0.612061168992 -0.280383142359 0.739436555016
-0.887250003871 0.258893773768 0.38178717178
-0.0971745158475 -0.994342015264 -0.0429076933695
-0.120898756509 -0.759794654167 -0.638823586113
-0.0208969821213 0.489420208746 -0.871797668848
0.482396765336 -0.834880934468 -0.265079584379
0.66383194755 0.726669941038 0.176855710123
from
>>>nfaces.py 32.txt 32fa.txt
here's a sample of the output
point vertex_0_2_12 -0.0288974102653 -0.851924110364 0.66948446135
point vertex_13_24_27 -0.122445373457 1.00045419254 0.398644871129
point vertex_0_7_15 -0.482610585141 -0.963539328269 0.1161816686
point vertex_3_10_17 0.848541556491 -0.639691741968 -0.213520247975
point vertex_1_6_14 -1.05498774772 0.248634080415 -0.00106786881656
point vertex_13_24_31 -0.291794484064 0.826881428947 0.637135994637
point vertex_1_6_25 -1.07023716261 -0.0479998647263 0.164643927545
point vertex_3_17_28 1.05498774772 -0.248634080415 0.00106786881656
point vertex_1_15_25 -0.975134617659 -0.4675255693 -0.073154040374
point vertex_4_8_10 0.291794484064 -0.826881428947 -0.637135994637
point vertex_9_19_20 -0.400242088946 0.638705226984 -0.778897359983
normal face_30 -0.617395768043 -0.359222235879 -0.699844161834
face face_31 vertex_21_29_31
face face_31 vertex_13_21_31
face face_31 vertex_13_24_31
face face_31 vertex_6_24_31
face face_31 vertex_6_29_31
normal face_31 -0.426352071227 0.457758408728 0.780180203927
As you can see, there are exactly 12 significant figures.
BTW, I did not write either script, I can post them, if you think
doing so may show where it may be possible to escape to "decimal"
as you say, perhaps?
"Decimal" has a specific meaning in Python different to how you
appear to be using it. It looks to me like you are working with
floats rather than decimals.
Right.
from math import pi, asin, atan2, cos, sin, sqrt
from math import pi, asin, atan2, cos, sin, sqrt
Two different scripts.
"crosspoint" appears to be a custom Python program that we know
nothing about, apart from what you tell us.
I posted it in it's entirety.
"==" is used for equality testing:
x == 2
returns False, because x does not equal 2; but
x == 1
returns True, because x does currently equal 1.
Ah! Thanks!
Almost. The "from math import blah..." line extracts a bunch of maths
functions from the math library and makes them available to your
program.
I'm seeing this now. Python does minimal things, extras are called in
only as needed.
The number of significant figures is more fundamental that that; your
operating system understands about floating point numbers (so-called
"C doubles") with 53 *binary* significant figures, or about 17 *decimal*
figures. So I'm not sure why you think there are only 12.
Sure. But see the above, from Py 2.7.2, it's precisely 12 sig.figs.!
You need to show us how the output is generated in the first place.
Done, above.
Not a chance without some major changes to your program.
I was afraid of this.
I can't say how major without seeing more of your program.
Well, OK, here's the first, mkpoints.py
#!/usr/bin/env python
# Distribute a set of points randomly across a sphere, allow them
# to mutually repel and find equilibrium.
import sys
import string
import random
from math import pi, asin, atan2, cos, sin, sqrt
args = sys.argv[1:]
if len(args) > 0:
n = string.atoi(sys.argv[1])
args = args[1:]
else:
n = 7
if len(args) > 0:
outfile = open(args[0], "w")
args = args[1:]
else:
outfile = sys.stdout
points = []
def realprint(a):
for i in range(len(a)):
outfile.write(str(a[i]))
if i < len(a)-1:
outfile.write(" ")
else:
outfile.write("\n")
def pprint(*a):
realprint(a)
for i in range(n):
# Invent a randomly distributed point.
#
# To distribute points uniformly over a spherical surface, the
# easiest thing is to invent its location in polar coordinates.
# Obviously theta (longitude) must be chosen uniformly from
# [0,2*pi]; phi (latitude) must be chosen in such a way that
# the probability of falling within any given band of latitudes
# must be proportional to the total surface area within that
# band. In other words, the probability _density_ function at
# any value of phi must be proportional to the circumference of
# the circle around the sphere at that latitude. This in turn
# is proportional to the radius out from the sphere at that
# latitude, i.e. cos(phi). Hence the cumulative probability
# should be proportional to the integral of that, i.e. sin(phi)
# - and since we know the cumulative probability needs to be
# zero at -pi/2 and 1 at +pi/2, this tells us it has to be
# (1+sin(phi))/2.
#
# Given an arbitrary cumulative probability function, we can
# select a number from the represented probability distribution
# by taking a uniform number in [0,1] and applying the inverse
# of the function. In this case, this means we take a number X
# in [0,1], scale and translate it to obtain 2X-1, and take the
# inverse sine. Conveniently, asin() does the Right Thing in
# that it maps [-1,+1] into [-pi/2,pi/2].
theta = random.random() * 2*pi
phi = asin(random.random() * 2 - 1)
points.append((cos(theta)*cos(phi), sin(theta)*cos(phi), sin(phi)))
# For the moment, my repulsion function will be simple
# inverse-square, followed by a normalisation step in which we pull
# each point back to the surface of the sphere.
while 1:
# Determine the total force acting on each point.
forces = []
for i in range(len(points)):
p = points[i]
f = (0,0,0)
ftotal = 0
for j in range(len(points)):
if j == i: continue
q = points[j]
# Find the distance vector, and its length.
dv = (p[0]-q[0], p[1]-q[1], p[2]-q[2])
dl = sqrt(dv[0]**2 + dv[1]**2 + dv[2]**2)
# The force vector is dv divided by dl^3. (We divide by
# dl once to make dv a unit vector, then by dl^2 to
# make its length correspond to the force.)
dl3 = dl ** 3
fv = (dv[0]/dl3, dv[1]/dl3, dv[2]/dl3)
# Add to the total force on the point p.
f = (f[0]+fv[0], f[1]+fv[1], f[2]+fv[2])
# Stick this in the forces array.
forces.append(f)
# Add to the running sum of the total forces/distances.
ftotal = ftotal + sqrt(f[0]**2 + f[1]**2 + f[2]**2)
# Scale the forces to ensure the points do not move too far in
# one go. Otherwise there will be chaotic jumping around and
# never any convergence.
if ftotal > 0.25:
fscale = 0.25 / ftotal
else:
fscale = 1
# Move each point, and normalise. While we do this, also track
# the distance each point ends up moving.
dist = 0
for i in range(len(points)):
p = points[i]
f = forces[i]
p2 = (p[0] + f[0]*fscale, p[1] + f[1]*fscale, p[2] + f[2]*fscale)
pl = sqrt(p2[0]**2 + p2[1]**2 + p2[2]**2)
p2 = (p2[0] / pl, p2[1] / pl, p2[2] / pl)
dv = (p[0]-p2[0], p[1]-p2[1], p[2]-p2[2])
dl = sqrt(dv[0]**2 + dv[1]**2 + dv[2]**2)
dist = dist + dl
points[i] = p2
# Done. Check for convergence and finish.
sys.stderr.write(str(dist) + "\n")
if dist < 2.10005e-16:
break
# Output the points.
for x, y, z in points:
pprint(x, y, z)
The nfaces.py used on the output from this is a tad longer!
Well, yes, but only with some significant changes
Well, it seems I may have to learn "decimal" python, and insert it
into the above, and nfaces.py.
John.
_______________________________________________
Tutor maillist - Tutor@python.org
To unsubscribe or change subscription options:
http://mail.python.org/mailman/listinfo/tutor
