Geodesic curve on a torus

Requires Rhino5, GH, and GhPython.

```import rhinoscriptsyntax as rs
import Rhino.Geometry as rg
import random
import math
def residual(_x):
x=_x.X
y=_x.Y
z=_x.Z
L=math.sqrt(x**2+z**2)
return ((x-R*x/L)**2)+((z-R*z/L)**2)+(y**2)-r**2
def jacob(_x):
x=_x.X
y=_x.Y
z=_x.Z
L=math.sqrt(x**2+z**2)
fx=2*(x-R*x/L)*(1-R/L+R*x/(L**3)*x)+2*(z-R*z/L)*(R*z/(L**3)*x)
fz=2*(z-R*z/L)*(1-R/L+R*z/(L**3)*z)+2*(x-R*x/L)*(R*x/(L**3)*z)
fy=2*y
return rg.Vector3d(fx,fy,fz)
def pinv(_x):
x=_x.X
y=_x.Y
z=_x.Z
N=x**2+y**2+z**2
return rg.Vector3d(x/N,y/N,z/N)
def varphi(_x):
J=jacob(X)
iJ=pinv(J)
N=_x.Length
lam=_x*J
_x=_x-lam*iJ
_x.Unitize()
_x*=N
return _x
def psi(_x):
for i in range(10):
J=jacob(_x)
iJ=pinv(J)
Resid=residual(_x)
dx=iJ*Resid
_x-=dx
return _x
dt=0.1
R=10
r=3
plane=rs.WorldZXPlane()
b=rg.Torus(plane,R,r)
if x:
for t in range(10):
V=varphi(V)
X+=V*dt
X=psi(X)
g=f
else:
X=rg.Point3d(13,0.1,0)
V=rg.Vector3d(-0.15,0.9,0.3)
f=rg.Polyline()
```

page revision: 5, last edited: 08 Dec 2013 12:49