Author: Sachi Hashimoto
Views : 137
from sage.ext.fast_eval import fast_float

var('x,y,z,a,b,c')
f = y^4 + 5*x^4 - 6*x^2*y^2 + 6*x^3*z + 26*x^2*y*z + 10*x*y^2*z - 10*y^3*z - 32*x^2*z^2 - 40*x*y*z^2 + 24*y^2*z^2 + 32*x*z^3 - 16*y*z^3
ff = fast_float(f,'x','y','z')

gx_ff = fast_float(gradx, 'x', 'y', 'z')
gy_ff = fast_float(grady, 'x', 'y', 'z')
gz_ff = fast_float(gradz, 'x', 'y', 'z')

pos_sphere = sqrt(16-x^2-y^2)
pos_ff = fast_float(pos_sphere,'x','y')

(x, y, z, a, b, c)

%time
N = 1000
R = 4
TOL = 1

# Carve up the xy plane into a (2N+1) x (2N+1) grid and try to find a z value of an intersectin at each lattice point
our_points = []
for i in range(-N,N+1):
if not i%100:
print i
x0 = R*float(i)/N
x_used = False
recent_y = None
for j in range(-N,N+1):
y0 = R*float(j)/N
# Compute the z value on the sphere
pos_z = pos_ff(x0,y0)
# Compute the function value and the gradient.
# If the gradient suggests we can get to an actual solution without changing by more than the size of our lattice, do it.
pos_val = ff(x0,y0,pos_z)
if abs(N*pos_ratio) < R/2:
# Since we're not too far off, adjust our point to lie on our surface better and consider it done.
# Negatives now!
neg_z = -pos_z
neg_val = ff(x0,y0,neg_z)
if abs(N*neg_ratio) < R/2:

print len(our_points)

-1000 -900 -800 -700 -600 -500 -400 -300 -200 -100 0 100 200 300 400 500 600 700 800 900 1000 7102 CPU time: 131.62 s, Wall time: 139.19 s
euclidean_distance = fast_float(sqrt((x-a)^2+(y-b)^2+(z-c)^2),'x','y','z','a','b','c') # There's probably something built in, but this'll do

# Now we attempt to put the points in neighborly order

MAX_DIST = 2
points_copy = our_points[:]
neighbors = [[points_copy.pop(3000)]] # Start with not the first point in case the corner is weird
counter = 0
while points_copy:
if not counter%1000:
print counter
counter += 1
cur_point = neighbors[-1][-1]
# Find the nearest point of the remaining ones
nn = min(range(len(points_copy)), key=lambda i: euclidean_distance(*(points_copy[i]+cur_point)))
dist = euclidean_distance(*(points_copy[nn]+cur_point))
if dist < MAX_DIST:
# If it's not out of range, let's use it
neighbors[-1].append(points_copy.pop(nn))
else:
# Otherwise, check if it should loop...
if euclidean_distance(*(neighbors[-1][-1] + neighbors[-1])) < MAX_DIST:
neighbors[-1].append(neighbors[-1])
# ...and regardless start a new chain
neighbors.append([points_copy.pop(nn)])
# Final loop check
if euclidean_distance(*(neighbors[-1][-1] + neighbors[-1])) < MAX_DIST:
neighbors[-1].append(neighbors[-1])

print len(neighbors), [len(block) for block in neighbors]

0 1000 2000 3000 4000 5000 6000 7000 2 [3552, 3552]

A = implicit_plot3d(f==0, (x,-4,4), (y,-4,4), (z,-4,4), plot_points=100, frame=false, opacity=.5) #Plot the surface
C = implicit_plot3d(x^2+y^2+z^2==16, (x,-4,4), (y,-4,4), (z,-4,4), plot_points=80, frame=false, color="seashell", opacity=.2)

the_visualization= sum((line(block, color="black") for block in neighbors), A+C) #Plot blocks of neighbors

show(the_visualization.rotate((-1,-1,3),-5*pi/12).scale(2))


Error in lines 1-1 Traceback (most recent call last): File "/cocalc/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 1188, in execute flags=compile_flags) in namespace, locals File "", line 1, in <module> NameError: name 'f' is not defined