Masaaki MIKI

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()