Curvature Calculator

var('t')                 ##declare the variables you'll use below.
r(t)=[(1*t),(1*t)^2]     ##Give a parametrization of the curve
tBounds=(t,-2,2)         ##State the bounds you want in your plot
point=1                  ##Give a t-value at which you would like to draw all vectors.

##Computation time
v=diff(r,t)              #compute the velocity
a=diff(v,t)              #compute the acceleration
speed = v.norm()         #.norm() computes the length of a vector.
T = v/speed              #compute the unit tangent vector
dTdt = diff(T,t)         #compute dT/dt
dTds = diff(T,t)/speed   #compute dT/ds (the curvatuve vector)
curvature = dTds.norm()  #the curvatuve is the length of the curvatuve vector.
N = dTdt/dTdt.norm()     #The unit normal vector

##Print the results.  I added .simplify_full().simplify_trig() to the end of each, so that the results will look simplified.
table([
[r"$\vec r$",r(t).simplify_full().simplify_trig()],
[r"$\vec v$",v(t).simplify_full().simplify_trig()],
[r"$\vec a$",a(t).simplify_full().simplify_trig()],
[r"$\frac{ds}{dt} = |\vec v| = $ speed",speed(t).simplify_full().simplify_trig()],
[r"$\vec T$",T(t).simplify_full().simplify_trig()],
[r"$\frac{d\vec T}{dt}$",dTdt(t).simplify_full().simplify_trig()],
[r"$\frac{d\vec T}{ds}$",dTds(t).simplify_full().simplify_trig()],
[r"$\left|\frac{d\vec T}{ds}\right| = \kappa $",curvature(t).simplify_full().simplify_trig()],
[r"$\vec N$",N(t).simplify_full().simplify_trig()],
])


##Print the results.  I added .simplify_full().simplify_trig() to the end of each, so that the results will look simplified.
print("At the time t = ",point, " we have the following:")
table([
[r"$\vec r$",r(point).simplify_full().simplify_trig()],
[r"$\vec v$",v(point).simplify_full().simplify_trig()],
[r"$\vec a$",a(point).simplify_full().simplify_trig()],
[r"$\frac{ds}{dt} = |\vec v| = $ speed",speed(point).simplify_full().simplify_trig()],
[r"$\vec T$",T(point).simplify_full().simplify_trig()],
[r"$\frac{d\vec T}{dt}$",dTdt(point).simplify_full().simplify_trig()],
[r"$\frac{d\vec T}{ds}$",dTds(point).simplify_full().simplify_trig()],
[r"$\left|\frac{d\vec T}{ds}\right| = \kappa $",curvature(point).simplify_full().simplify_trig()],
[r"$\vec N$",N(point).simplify_full().simplify_trig()],
])

##Now let's make a graph
p=parametric_plot(r,tBounds)                                        #this draws the curve
p+=plot(vector(r(point)), width=4, color='gray')
p+=plot(vector(v(point)), start=r(point), width=4, color='red')
p+=plot(vector(a(point)), start=r(point), width=4, color='green')
p+=plot(vector(T(point).simplify_full().simplify_trig()), start=r(point), width=4, color='black')
p+=plot(vector(N(point).simplify_full().simplify_trig()), start=r(point), width=4, color='purple')
p+=plot(vector(dTdt(point).simplify_full().simplify_trig()), start=r(point), width=4, color='yellow')
p+=plot(vector(dTds(point).simplify_full().simplify_trig()), start=r(point), width=4, color='blue')

print("Here's a plot of r, with the vectors v and a drawn at a point.")
print("gray-dashed = r, red = v, green = a, black = T, yellow = dT/dt, blue = dT/ds, purple=N")

show(p)

radiusOfCurvature=1/curvature(point)                    #compute the radius of curvature
centerOfCurvature=r(point)+radiusOfCurvature*N(point)   #compute the center of curvature

table([
[r"$\rho$",radiusOfCurvature.simplify_full().simplify_trig()],
[r"center of curvature", centerOfCurvature.simplify_full().simplify_trig()],
])

#Draw the circle of curvature on top of everything else.
p+=parametric_plot(
    vector(centerOfCurvature.simplify_full().simplify_trig())
    +radiusOfCurvature.simplify_full().simplify_trig()
    *( vector(T(point).simplify_full().simplify_trig())*cos(t)
      +vector(N(point).simplify_full().simplify_trig())*sin(t))
    ,(t,0,2*pi),color='gray')
show(p)

print("If you are in 3D, then this will compute B and dB/ds for you. Otherwise, you'll get a lot of errors.")

B=T.cross_product(N)
dBds = B.diff(t)/speed

table([
[r"$N(t)$",B(t).simplify_full().simplify_trig()],
[r"$B(t)$",B(t).simplify_full().simplify_trig()],
[r"$\frac{d\vec B}{ds}(t)$", dBds(t).simplify_full().simplify_trig()],
[r"$\left|\frac{d\vec B}{ds}(t)\right|$", dBds(t).norm().simplify_full().simplify_trig()],
["",r"Remember to compare $\vec N$ and $d\vec B/ds$ to determine torsion"],
[r"$N(point)$",N(point).simplify_full().simplify_trig()],
[r"$B(point)$",B(point).simplify_full().simplify_trig()],
[r"$\frac{d\vec B}{ds}(point)$", dBds(point).simplify_full().simplify_trig()],
[r"$\left|\frac{d\vec B}{ds}(point)\right|$", dBds(point).norm().simplify_full().simplify_trig()],
])