CoCalc Public FilesMotion in 3D.sagews
Author: Thomas Krainer
Views : 27
Description: Motion in 3-dimensional space
# Motion in 3D
# Author: Thomas Krainer (Penn State Altoona)
# Version: 02/11/2016
# Prior Version: 08/26/2011
#
# Changes: Norms changed from method .norm() to calculations by hand to fix
# issues with running the program.
#
# Based on "3-D motion and vector calculus" by Robert A. Beezer
# See http://wiki.sagemath.org/interact/calculus
#

t=var('t')
u,v=var('u,v')

output='<h1>Motion in 3-dimensional space</h1>'

output+='\\begin{aligned} &\\textrm{Position:} &&\\vec{p}(t) = \\langle x(t),y(t),z(t) \\rangle' output+=' \\quad \\textrm{for } t \in [t_{\\min},t_{\\max}] \\\ ' output+='&\\textrm{Velocity:} &&\\vec{v} = \\frac{d\\vec{p}}{dt} \\\ ' output+='&\\textrm{Speed:} &&v = |\\vec{v}| \\\ ' output+='&\\textrm{Acceleration:} &&\\vec{a} = \\frac{d\\vec{v}}{dt} = \\frac{d^2\\vec{p}}{dt^2} = a_{T}\\vec{T} + a_N\\vec{N} = \\left(\\frac{dv}{dt}\\right) \\vec{T} + \\left(\\kappa v^2\\right) \\vec{N} = \\left(\\frac{\\vec{v}\\cdot\\vec{a}}{v}\\right)\\vec{T} + \\left(\\frac{|\\vec{v}\\times\\vec{a}|}{v}\\right)\\vec{N} \\\ ' output+='&\\textrm{Unit tangent:} &&\\vec{T} = \\frac{\\vec{v}}{v} \\\ ' output+='&\\textrm{Unit normal:} &&\\vec{N} = \\frac{\\frac{d\\vec{T}}{dt}}{\\left|\\frac{d\\vec{T}}{dt}\\right|} \\\ ' output+='&\\textrm{Binormal:} &&\\vec{B} = \\vec{T}\\times\\vec{N} \\\ ' output+='&\\textrm{Curvature:} &&\\kappa = \\left|\\frac{d\\vec{T}}{ds}\\right| = \\frac{\\left|\\frac{d\\vec{T}}{dt}\\right|}{v} = \\frac{|\\vec{v}\\times\\vec{a}|}{v^3} \\\ ' output+='&\\textrm{Torsion:} &&\\tau = \\frac{(\\vec{v}\\times\\vec{a})\\cdot \\frac{d^3\\vec{p}}{dt^3}}{|\\vec{v}\\times\\vec{a}|^2} \\\ ' output+='&\\textrm{Frenet-Serret formulas:} &&\\begin{cases} \\frac{d\\vec{T}}{ds} = \\kappa\\vec{N} \\\ \\frac{d\\vec{N}}{ds} = -\\kappa\\vec{T} + \\tau\\vec{B} \\\ \\frac{d\\vec{B}}{ds} = -\\tau\\vec{N} \\end{cases} \\\ ' output+='&\\textrm{Osculating circle:} &&\\textrm{Circle of radius}\\; 1/\\kappa \\; \\textrm{in the plane through}\\;\\vec{p}\\;\\textrm{with normal vector}\\;\\vec{B}\\textrm{, centered at}\\; \\vec{p} + \\bigl(1/\\kappa\\bigr)\\vec{N}. ' output+=' \\end{aligned}'

html(output)

@interact
def Motion(curve=input_box(default='[(1/15)*(t^2+1)*cos(t),(1/15)*(t^2+1)*sin(t),t]',type=str,width=100,label='$[x(t),y(t),z(t)]=$'),\
interval=input_box(default='[0,3*pi]',type=str,width=30,label='$[t_{\\min},t_{\\max}]=$'),\
param=slider(0,100,1,label='Parameter $t_0=$',display_value=False),\
show_position=('Show position (blue point):',True),\
show_velocity=('Show velocity (green):',False),\
show_acceleration=('Show acceleration (pink) with tangential (green) and normal (red) projections:',False),\
show_TNBframe=('Show TNB-frame (green/red/blue):',False),\
show_osccircle=('Show osculating circle:',False),\
show_computations=('Show computations:',True)):

[tmin,tmax]=sage_eval(interval)
ndigits=5 # Numerical precision up to 5 digits behind the decimal dot

#
# Vector and scalar functions pertaining to the motion
#
position(t)=sage_eval(curve,locals={'t':t})
velocity = derivative(position,t)
acceleration = derivative(velocity,t)
jerk = derivative(acceleration,t)
speed = sqrt(velocity[0]^2+velocity[1]^2+velocity[2]^2)
tangent = (1/speed)*velocity
dT = derivative(tangent,t)
dTnorm = sqrt(dT[0]^2+dT[1]^2+dT[2]^2)
normal = (1/(dTnorm))*dT
binormal = tangent.cross_product(normal)
vcrossa=velocity.cross_product(acceleration)
vcrossanorm=sqrt(vcrossa[0]^2+vcrossa[1]^2+vcrossa[2]^2)
curvaturefunction = (1/speed^3)*vcrossanorm
torsionfunction = (1/(vcrossanorm^2))*(vcrossa.dot_product(jerk))

picture=parametric_plot3d(position(t),(t,tmin,tmax),color='red',thickness=3)
plotsize=(vector(picture.bounding_box()[1])-vector(picture.bounding_box()[0])).norm()

#
# Scalar quantities at point
#
t0=RR(tmin+(param/100)*(tmax-tmin))
pos_tzero = position(t0)
tangent_comp = (1/speed(t0))*velocity(t0).dot_product(acceleration(t0))
normal_comp = (1/speed(t0))*vcrossanorm(t0)
curvature = curvaturefunction(t0)
torsion = torsionfunction(t0)

#
# Vector, plane, and circle plot
#
pos = point3d(pos_tzero, color='blue',size=10)
vel = arrow3d(pos_tzero, pos_tzero + velocity(t0), rgbcolor=(0,0.5,0))
tan = arrow3d(pos_tzero, pos_tzero + tangent(t0), rgbcolor=(0,1,0) )
nor = arrow3d(pos_tzero, pos_tzero + normal(t0), rgbcolor=(0.5,0,0))
bin = arrow3d(pos_tzero, pos_tzero + binormal(t0), rgbcolor=(0,0,0.5))
acc = arrow3d(pos_tzero, pos_tzero + acceleration(t0), rgbcolor=(1,0,1))
tanproj = arrow3d(pos_tzero, pos_tzero + tangent_comp*tangent(t0), rgbcolor=(0,1,0))
norproj = arrow3d(pos_tzero, pos_tzero + normal_comp*normal(t0), rgbcolor=(0.5,0,0))
TNplane = parametric_plot3d(pos_tzero+u*tangent_comp*tangent(t0)+v*normal_comp*normal(t0),\
(u,0,1),(v,0,1),opacity=0.2,color='yellow')
oscdisk = parametric_plot3d(pos_tzero+(1/curvature)*(1+v*sin(u))*normal(t0)\
+(1/curvature)*v*cos(u)*tangent(t0),(u,0,2*pi),(v,0,1),color='lightblue',thickness=3,opacity=0.5)
osccircle = parametric_plot3d(pos_tzero+(1/curvature)*(1+sin(u))*normal(t0)\
+(1/curvature)*cos(u)*tangent(t0),(u,0,2*pi),color='darkblue',thickness=3,opacity=0.8)

#
# Assemble plot
#
if (show_position == True):
picture+=pos
if (show_velocity == True):
picture+=vel
if (show_acceleration == True):
picture+=acc+tanproj+norproj+TNplane
if (show_TNBframe == True):
picture+=tan+nor+bin
if (show_osccircle == True):
picture+=oscdisk+osccircle

#
# Assemble written output
#
output='\\begin{aligned} ' output+='&\\textrm{Position function:} &&\\vec{p}(t) = \\left\\langle %s,%s,%s \\right\\rangle \\quad \\textrm{for}\\; t \\in [%s,%s] \\\ '%(latex(position[0](t)),latex(position[1](t)),latex(position[2](t)),latex(tmin),latex(tmax)) output+='&\\textrm{Parameter value:} &&t_0 = %s \\\ '%(latex(t0.n(digits=ndigits))) output+='&\\textrm{Position:} &&\\vec{p}(t_0) = \\left\\langle %s,%s,%s \\right\\rangle \\\ '%(latex(pos_tzero[0].n(digits=ndigits)),latex(pos_tzero[1].n(digits=ndigits)),latex(pos_tzero[2].n(digits=ndigits))) output+='&\\textrm{Velocity:} &&\\vec{v}(t_0) = \\left\\langle %s,%s,%s \\right\\rangle \\\ '%(latex(velocity[0](t0).n(digits=ndigits)),latex(velocity[1](t0).n(digits=ndigits)),latex(velocity[2](t0).n(digits=ndigits))) output+='&\\textrm{Speed:} &&v(t_0) = %s \\\ '%(latex(speed(t0).n(digits=ndigits))) output+='&\\textrm{Acceleration:} &&\\vec{a}(t_0) = \\left\\langle %s,%s,%s \\right\\rangle = %s \\cdot\\vec{T}(t_0) + %s \\cdot\\vec{N}(t_0) \\\ '%(latex(acceleration[0](t0).n(digits=ndigits)),latex(acceleration[1](t0).n(digits=ndigits)),latex(acceleration[2](t0).n(digits=ndigits)),latex(tangent_comp.n(digits=ndigits)),latex(normal_comp.n(digits=ndigits))) output+='&\\textrm{Unit tangent:} &&\\vec{T}(t_0) = \\left\\langle %s,%s,%s \\right\\rangle \\\ '%(latex(tangent[0](t0).n(digits=ndigits)),latex(tangent[1](t0).n(digits=ndigits)),latex(tangent[2](t0).n(digits=ndigits))) output+='&\\textrm{Unit normal:} &&\\vec{N}(t_0) = \\left\\langle %s,%s,%s \\right\\rangle \\\ '%(latex(normal[0](t0).n(digits=ndigits)),latex(normal[1](t0).n(digits=ndigits)),latex(normal[2](t0).n(digits=ndigits))) output+='&\\textrm{Binormal:} &&\\vec{B}(t_0) = \\left\\langle %s,%s,%s \\right\\rangle \\\ '%(latex(binormal[0](t0).n(digits=ndigits)),latex(binormal[1](t0).n(digits=ndigits)),latex(binormal[2](t0).n(digits=ndigits))) output+='&\\textrm{Curvature:} &&\\kappa(t_0) = %s \\\ '%(latex(curvature.n(digits=ndigits))) output+='&\\textrm{Torsion:} &&\\tau(t_0) = %s '%(latex(torsion.n(digits=ndigits))) output+=' \\end{aligned}'

#
# Show written output and picture
#
if (show_computations == True):
html(output)
show(picture,aspect_ratio=[1,1,1])


# Motion in 3-dimensional space

\begin{aligned} &\textrm{Position:} &&\vec{p}(t) = \langle x(t),y(t),z(t) \rangle \quad \textrm{for } t \in [t_{\min},t_{\max}] \\ &\textrm{Velocity:} &&\vec{v} = \frac{d\vec{p}}{dt} \\ &\textrm{Speed:} &&v = |\vec{v}| \\ &\textrm{Acceleration:} &&\vec{a} = \frac{d\vec{v}}{dt} = \frac{d^2\vec{p}}{dt^2} = a_{T}\vec{T} + a_N\vec{N} = \left(\frac{dv}{dt}\right) \vec{T} + \left(\kappa v^2\right) \vec{N} = \left(\frac{\vec{v}\cdot\vec{a}}{v}\right)\vec{T} + \left(\frac{|\vec{v}\times\vec{a}|}{v}\right)\vec{N} \\ &\textrm{Unit tangent:} &&\vec{T} = \frac{\vec{v}}{v} \\ &\textrm{Unit normal:} &&\vec{N} = \frac{\frac{d\vec{T}}{dt}}{\left|\frac{d\vec{T}}{dt}\right|} \\ &\textrm{Binormal:} &&\vec{B} = \vec{T}\times\vec{N} \\ &\textrm{Curvature:} &&\kappa = \left|\frac{d\vec{T}}{ds}\right| = \frac{\left|\frac{d\vec{T}}{dt}\right|}{v} = \frac{|\vec{v}\times\vec{a}|}{v^3} \\ &\textrm{Torsion:} &&\tau = \frac{(\vec{v}\times\vec{a})\cdot \frac{d^3\vec{p}}{dt^3}}{|\vec{v}\times\vec{a}|^2} \\ &\textrm{Frenet-Serret formulas:} &&\begin{cases} \frac{d\vec{T}}{ds} = \kappa\vec{N} \\ \frac{d\vec{N}}{ds} = -\kappa\vec{T} + \tau\vec{B} \\ \frac{d\vec{B}}{ds} = -\tau\vec{N} \end{cases} \\ &\textrm{Osculating circle:} &&\textrm{Circle of radius}\; 1/\kappa \; \textrm{in the plane through}\;\vec{p}\;\textrm{with normal vector}\;\vec{B}\textrm{, centered at}\; \vec{p} + \bigl(1/\kappa\bigr)\vec{N}. \end{aligned}