# Set x, y, z as functions of u parametrizing your space curve. u,v = var('u,v'); # Formula for unknot: x = cos(u); y = sin(u); z = 0; # Formulas for p,q torus knot: # If p=2 and q=3, then this is the trefoil knot. #p=2; q=3; #r = cos(q*u)+2; #x = r*cos(p*u); #y = r*sin(p*u); #z = -sin(q*u); # Formula for figure eight: #x = (2+cos(2*u))*cos(3*u); #y = (2+cos(2*u))*sin(3*u); #z = sin(4*u); # gamma is now set to be a vector representation of your space curve. # gamma_T, gamma_N, gamma_B form a Frenet frame for gamma (assuming some non-vanishing conditions on derivatives of gamma). gamma = vector([x,y,z]); dgamma = gamma.diff(u) gamma_T = dgamma.normalized(); gamma_N = gamma_T.diff(u).normalized(); gamma_B = gamma_T.cross_product(gamma_N); # radius of tube around space curve radius = 0.2; # Set a, b so that the space curve is parametrized over [a,b]. # Set n to be the number of subdivisions of [0,2*pi] (for the tube). # Set k to be the number of subdivisions of [a,b]. a = 0; b = 2*pi; n = 6; def curve_len(gamma,a,b): return numerical_integral(gamma.diff(u).norm(),a,b)[0]; approx_k = (n*curve_len(gamma,a,b)/(sqrt(3)*pi*radius)).n(); def nearest_even(x): f = floor(x); if f%2 == 0: return f; else: return f+1; k = nearest_even(approx_k); # This k will produce approximately equilateral faces in the triangulation. # The closer gamma is to a unit speed curve, the better the approximation. ustep = (b-a)/k; vstep = 2*pi/n; # tube gives a Frenet tube around gamma # ribbon = normal_ribbon gives the normal ribbon on gamma # binormal_ribbon gives the binormal ribbon on gamma tube = gamma + radius*(cos(v)*gamma_N + sin(v)*gamma_B); normal_ribbon = gamma + v*gamma_N; binormal_ribbon = gamma + v*gamma_B; ribbon = normal_ribbon; # points is a list of lists of vectors; its i-th entry is a list of n vectors to be plotted. points = [[tube(ustep*i,vstep*(j+(i%2)/2)).n(32) for j in [0..n-1]] for i in [0..k]]; def tri_normal(v0,v1,v2): return (v1-v0).cross_product(v2-v0).normalized(); def write_face(v0,v1,v2): normal = tri_normal(v0,v1,v2); out = " facet normal " + repr(normal[0]) + " " + repr(normal[1]) + " " + repr(normal[2]) + "\n"; out += " outer loop\n" out += " vertex " + repr(v0[0]) + " " + repr(v0[1]) + " " + repr(v0[2]) + "\n" out += " vertex " + repr(v1[0]) + " " + repr(v1[1]) + " " + repr(v1[2]) + "\n" out += " vertex " + repr(v2[0]) + " " + repr(v2[1]) + " " + repr(v2[2]) + "\n" out += " endloop\n" out += " endfacet\n" return out; def list_to_stl(list, name): i_len = len(list); j_len = len(list[0]); out = "solid " + name + "\n"; for i in [0..i_len-2]: if i%2 == 0: for j in [0..j_len-1]: v0 = list[i][j]; v1 = list[i][(j+1)%j_len]; v2 = list[i+1][j]; out += write_face(v0,v1,v2); v0 = list[i+1][j]; v1 = list[i][(j+1)%j_len]; v2 = list[i+1][(j+1)%j_len]; out += write_face(v0,v1,v2); else: for j in [0..j_len-1]: v0 = list[i][j]; v1 = list[i][(j+1)%j_len]; v2 = list[i+1][(j+1)%j_len]; out += write_face(v0,v1,v2); v0 = list[i][j]; v1 = list[i+1][(j+1)%j_len]; v2 = list[i+1][j]; out += write_face(v0,v1,v2); out += "endsolid " + name + "\n"; return out; f = file('/Users/ormsbyk/Dropbox/math/3d/space_curve/unknot.stl','w'); # Replace the first string with your desired file location / name. f.write(list_to_stl(points, "unknot")); # Replace "unknot" with a good name for your knot / curve. f.close();