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