I'm working on visualizing a flow in R^3, and specifically focusing on how the flow behaves near a torus on inner radius 1 and outer radius 2, centered at the origin. The idea is to take a point (x,y,z), and look at this point using coordinates (r,z), in the plane containing (x,y,z) and the origin. Then I take equations in polar coordinates in this plane, and apply them to the point. Moving back into R^3, I take a rotation of the point, and add in the result of the flow restricted to the plane. Most of these conditions also require that I know the minimum distance from (x,y,z) to the torus mentioned.
This is taking a very long time, especially when I try to animate the flow. The last animation I tried took 20 hours, and that was before I found a flaw in the construction. This one works, but I'm hoping to optimize it before I try it again. Any suggestions on how to speed things up would be appreciated.
x,y,u,v = var('x y u v')
ep = 0.1 #step size
length = 350 #iterations
#initial conditions
omega_seeds = [(-1,-1,2),(1,-1,2),(-1,1,2),(1,1,2),(2,2,2),(-2,2,2),(2,-2,2),(-2,-2,2)]
twopi = 2 * RDF.pi()
c(r) = 2/3*(4**(r+1)-4)
def norm2((x,y)):
"""
Norm of `(x,y)`
"""
return float(sqrt((x^2 + y^2)))
def norm3((x,y,z)):
"""
Norm squared of `(x,y,z)`
"""
return x^2 + y^2 + z^2
def norm((x,y,z)):
"""
Norm of `(x,y,z)`
"""
return float(sqrt(norm3((x,y,z))))
def dist((x,y,z),(xx,yy,zz)):
"""
Distance from `(x,y,z)` to `(xx,yy,zz)`
"""
return norm((x - xx, y - yy,z-zz))
def atan2((x,y)): #generalized arctangent
if x == 0 and y < 0:
return 3*pi/2
elif x == 0 and y > 0:
return pi/2
else:
a = float(arctan(abs(y/x)))
if x > 0:
if y > 0:
return a
elif y < 0:
return -a
elif y == 0:
return 0
elif x < 0:
if y > 0:
return pi-a
elif y < 0:
return a-pi
#Functions to determine distances of a point from a specified torus
def Torus_dist((x,y,z),n):
A = vector([0,c(n),0])
if n%2 == 0:
B = vector([x,y,0])
else:
B = vector([0,y,z])
if A == B:
return 4^n
else:
t = (2*4^n)/(norm(A - B))
v = t*B+ (1-t)*A
C = vector([v[0],v[1],v[2]])
if dist(C,(x,y,z)) <= 4^n:
return 0
else:
D = vector([x,y,z])
t2 = 4^n/norm(C-D)
v2 = t2*D + (1-t2)*C
E = vector([v2[0],v2[1],v2[2]])
return dist((x,y,z),E)
def dz(x,y): #vertical component of the flow restricted to the plane
r = norm2((x,y))
if r <=1:
return 0
if r >=2:
return y
theta = atan2((x,y))
dr = cos(theta)*(1-1/r^2)
dtheta = -sin(theta)*(1+1/r^2)
return float(dr*cos(theta)-r*sin(theta)*dtheta)
def dr(x,y): #horizontal component of the flow restricted to the plane
r = norm2((x,y))
if r <=1:
return 0
if r >=2:
return 0
theta = atan2((x,y))
dr = cos(theta)*(1-1/r^2)
dtheta = -sin(theta)*(1+1/r^2)
return float(dr*sin(theta) + r*cos(theta)*dtheta)
def W1((x,y,z)): #dx for the points in R^3
r = 2 - norm2((x,y))
T = max(0,1-Torus_dist((x,y,z),0))
theta = atan2((x,y))
if 0<=z<=2:
return T*y + dr(r,z)*cos(theta)
elif -2<=z<=0:
return (-1)*T*y + dr(r,z)*cos(theta)
else:
return 0
def W2((x,y,z)): #dy for the points in R^3
r = 2 - norm2((x,y))
T = max(0,1-Torus_dist((x,y,z),0))
theta = atan2((x,y))
if 0<=z<=2:
return (-1)*T*x + dr(r,z)*sin(theta)
elif -2<=z<=0:
return T*x + dr(r,z)*sin(theta)
else:
return 0
def W3((x,y,z)): #dz for the points in R^3
r = 2 - norm2((x,y))
if -2<=z<=2:
return -dz(r,z)
else:
return 0
def Plugplot((x,y,z)): #takes a point in R^3 and returns a point epsilon distance away, in the direction of the gradient at that point
n = norm((x - W1((x,y,z)), y - W2((x,y,z)),z-W3((x,y,z))))*(1/ep)
if n == 0:
return ((x,y,z))
else:
return ((x + W1((x,y,z))/n), (y + W2((x,y,z))/n), (z + W3((x,y,z))/n))
def PlugFlow((x,y,z)): #creates a collection of lines through the points generated by Plugplot
L = []
ic = (x,y,z)
if z >= -2:
for i in range(length):
L.append(line((ic,Plugplot(ic)),color='blue'))
ic = Plugplot(ic)
return L
"""
Draws the torus around which the flow exists
"""
a,b = 2,1 # outer and inner radii
Tx = (a + b*cos(u))*cos(v)
Ty = (a + b*cos(u))*sin(v)
Tz = b*sin(u)
T = parametric_plot3d([Tx,Ty,Tz], (u,0,twopi),(v,0,twopi), opacity=.2,aspect_ratio=1,color='black')
#Adds all of the to an array with the torus T
L = [T]
for i in omega_seeds:
L = L + PlugFlow(i)
show(sum(L) ,frame=true,aspect_ratio=1)
