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) f.Add(X) g=f a=rs.AddPolyline(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