In [1]:
import plotly.offline as py
from plotly.graph_objs import *
import numpy as np
from numpy import pi, cos, sin, exp, log, sqrt
py.init_notebook_mode()
In [2]:
def curve(rfun,tmin=-2,tmax=2,tpts=200,color='black'):
    domain = np.linspace(tmin,tmax,tpts)
    r = [[rfun(t)[i] for t in domain] for i in range(3)]
    trace = Scatter3d(x=r[0],y=r[1],z=r[2],mode='lines',
                      line=Line(color=color,width=3))
    return(trace)
In [3]:
def snake(u,v):
    lu, lv = len(u), len(v)
    path = []

    i, j = 0, 0
    istep, jstep = 1, 1
    while(i < lu):
        while(0 <= j < lv):
            path.append((u[i],v[j]))
            j += jstep
        j -= jstep
        i += istep
        jstep *= -1

    i -= istep
    i -= istep
    istep *= -1
    while(0 <= j < lv):
        while(0 <= i < lu):
            path.append((u[i],v[j]))
            i += istep
        i -= istep
        j += jstep
        istep *= -1
    return path

def surface(rfun,tmin=-2,tmax=2,tpts=20,umin=-2,umax=2,upts=20,color='black'):
    t, u = np.linspace(tmin,tmax,tpts), np.linspace(umin,umax,upts)
    path = snake(t,u)
    r = [[rfun(t,u)[i] for (t,u) in path] for i in (0,1,2)]
    trace = Scatter3d(x=r[0],y=r[1],z=r[2],mode='lines',
                      line=Line(color=color,width=3))
    return(trace)
In [4]:
def dot(u,v):
    '''3D dot product'''
    return u[0]*v[0]+u[1]*v[1]+u[2]*v[2]
def mag(v):
    '''magnitude of 3D vector'''
    return sqrt(dot(v,v))
def scale(c,v):
    '''scalar multiple c*v'''
    return [c*x for x in v]
def plus(u,v):
    return [u[i]+v[i] for i in (0,1,2)]
def minus(u,v):
    return [u[i]-v[i] for i in (0,1,2)]
def unit(v):
    '''direction (unit vector) of v'''
    return scale(1/mag(v),v)
def cross(u,v):
    '''3D cross product'''
    return [u[1]*v[2]-u[2]*v[1],u[2]*v[0]-u[0]*v[2],u[0]*v[1]-u[1]*v[0]]
In [5]:
dt=1e-6
def D(u,t0,dt=dt):
    '''approximate du/dt for u 3D vector'''
    u1 = u(t0+dt)
    u2 = u(t0-dt)
    du = minus(u2,u1)
    return scale(1/(dt+dt),du)
In [6]:
def T(r,t0,dt=dt):
    '''unit tangent vector'''
    rp = D(r,t0,dt=dt)
    return unit(rp)
def Tp(r,t0,dt=dt):
    '''dT/dt'''
    fun = lambda t: T(r,t,dt=dt)
    return D(fun, t0, dt=dt)
def kappa(r,t0,dt=dt):
    '''curvature'''
    mTp = mag(Tp(r,t0,dt=dt))
    mrp = mag(D(r,t0,dt=dt))
    return mTp/mrp
def N(r,t0,dt=dt):
    '''unit normal vector'''
    return unit(Tp(r,t0,dt=dt))
def B(r,t0,dt=dt):
    '''unit binormal vector'''
    return cross(T(r,t0,dt=dt),N(r,t0,dt=dt))
def Bp(r,t0,dt=dt):
    '''dB/dt'''
    fun = lambda t: B(r,t,dt=dt)
    return D(fun, t0, dt=dt)
def tau(r,t0,dt=dt):
    '''torsion'''
    mBp = mag(Bp(r,t0,dt=dt))
    mrp = mag(D(r,t0,dt=dt))
    return mBp/mrp
In [7]:
r = lambda t: [3*cos(2*t),3*sin(2*t),t]
helix=curve(r,color='red',tmin=0,tmax=4*pi)
py.iplot(Figure(data=Data([helix])))
In [8]:
t0=2
tangent=T(r,t0)
normal=N(r,t0)
binormal=B(r,t0)
curvature=kappa(r,t0)
radius=1/curvature
torsion=tau(r,t0)
print(tangent) #should be [-6*sin(4),6*cos(4),1] / sqrt(37)
print(normal) #should be [-cos(4),-sin(4),0]
print(binormal) #should be [2*cos(4),2*sin(4),0] / sqrt(37)
print(mag(tangent)) # should be 1
print(mag(normal)) # should be 1
print(mag(binormal)) #should be 1
print(dot(tangent,normal)) #should be 0
print(dot(normal,binormal)) #should be 0
print(dot(binormal,tangent)) #should be 0
print(curvature) #should be 12/37
print(radius) #should be 37/12
print(torsion) #should be 2/37
[-0.74650538292778879, 0.64475009595373911, -0.16439898730673402]
[0.65364390065944533, 0.75680225364739973, -2.2185647778222317e-06]
[0.12441609367123319, -0.10746005149818529, -0.98639392382941937]
1.0
1.0
1.0
2.67769691859e-11
-1.25606515748e-17
0.0
0.324323681881
3.08333944102
0.0344116740887

In [9]:
def osc_plane(curve,t0,dt=dt):
    T0=T(curve,t0,dt=dt)
    N0=N(curve,t0,dt=dt)
    base=curve(t0)
    arrow=lambda u,v: plus(scale(u,T0),scale(v,N0))
    return(lambda u,v: plus(base,arrow(u,v)))
def osc_circle(curve,t0,dt=dt):
    T0=T(curve,t0,dt=dt)
    N0=N(curve,t0,dt=dt)
    R=1/kappa(curve,t0,dt=dt)
    base=plus(curve(t0),scale(R,N0))
    arrow=lambda theta: plus(scale(cos(theta),T0),scale(sin(theta),N0))
    return(lambda theta: plus(base,scale(R,arrow(theta))))
In [11]:
arc=curve(r,color='red',tmin=0,tmax=4)
plane=surface(osc_plane(r,t0),color='yellow')
circ=curve(osc_circle(r,t0),color='blue',tmin=0,tmax=2*pi)
py.iplot(Figure(data=Data([arc,plane,circ])))
In []: