Initialization
| > | # First we must initialize...
restart; with(VectorCalculus): with(plots): # # osculatingCircles := proc(curve_eq,r,t_begin,t_end,numPts,nFrames,nLine,x_min,x_max,y_min,y_max) # osculatingCircles := proc(curve_eq,r,t_begin,t_end) local T,N,kappa,kappa_t0,cir_radius,cir_center,t_0,cc_x,cc_y,osc_circle, pt_x0,pt_y0,norm_line,graphX_min,graphY_min,graphX_max,graphY_max, graphX_length,graphY_length,t,opts,numPts,nFrames,nLine, x_min,x_max,y_min,y_max; opts:= [args[5..-1]]; # The normal line is shown by default. if nops(opts)<3 then nLine := true; else nLine := opts[3]; end if; # The default number of points (used by implicitplot) is 1000. if nops(opts)<1 then numPts := 1000; else numPts := opts[1]; end if; # The default number of frames shown is 50. if nops(opts)<2 then nFrames := 50; else nFrames := opts[2]; end if; # Here is the unit tangent vector T := simplify(diff(r(t),t)/Norm(diff(r(t),t))) assuming t::real; # Here is the unit normal vector N := simplify(diff(T,t)/Norm(diff(T,t))) assuming t::real; # We don't need the binormal. # Here is the curvature kappa := simplify(Curvature(r(t),t)) assuming t::real; # We can now find the equation of the osculating circle at r(t_0) # The osculating circle has radius 1/kappa # kappa_t0 is the curvature at t=t_0. # If the curvature is 0, then there is no osculating circle (it would # need to have an infinite radius). At such a point we should use a # tangent line. But instead, I'll use a circle with a really huge # radius (which looks the same when graphed). kappa_t0 := eval(kappa,t=t_0); cir_radius := piecewise(kappa_t0=0,999999999,1/kappa_t0); # We find its center by starting at r(t_0) on our curve and then # traveling (in the unit normal N(t_0) direction) 1/kappa # (the radius of the circle) units. cir_center := r(t_0) + cir_radius*eval(N,t=t_0); # Finally, we dot circle_center with e_x and e_y in order to # rip-off the x and y coordinates of the center. cc_x := cir_center . <1,0,0>; cc_y := cir_center . <0,1,0>; # Here is the equation of the osculating circle at r(t_0) osc_circle := (x-cc_x)^2 + (y-cc_y)^2 = cir_radius^2; # Here is a line from the curve (at r(t_0)) to the center of the circle. pt_x0 := r(t_0) . <1,0,0>; pt_y0 := r(t_0) . <0,1,0>; norm_line := simplify((pt_x0 - cc_x)*(y - pt_y0) = (pt_y0 - cc_y)*(x - pt_x0)); # Now let's plot some of these circles, normal lines, and the curve itself # To make sure that all of the circles fit in the plot window, we sample # some of the circles' extreme x,y coordinates (taking max's/min's). But # we also cap the length in case of crazy huge circles (like what happens # if we approach zero curvature) -- this is what the outer max's/min's and # graph??? variables are all about. graphX_min := min(seq(r(t_0).<1,0,0>,t_0=t_begin..t_end,(t_end-t_begin)/10)): graphY_min := min(seq(r(t_0).<0,1,0>,t_0=t_begin..t_end,(t_end-t_begin)/10)): graphX_max := max(seq(r(t_0).<1,0,0>,t_0=t_begin..t_end,(t_end-t_begin)/10)): graphY_max := max(seq(r(t_0).<0,1,0>,t_0=t_begin..t_end,(t_end-t_begin)/10)): graphX_length := graphX_max-graphX_min: graphY_length := graphY_max-graphY_min: if nops(opts)<7 then x_max := graphX_max+graphX_length*10: x_min := graphX_min-graphX_length*10: y_max := graphY_max+graphY_length*10: y_min := graphY_min-graphY_length*10: else x_min := opts[4]; x_max := opts[5]; y_min := opts[6]; y_max := opts[7]; end if; # to view this animation, first click on the plot, then some buttons will appear # where the "[Text][Math][C Maple Input]..." toolbar used to be (above) # # for smoother animations, set frames to be something like 100 # ***This calculation may take a minute or two*** if nLine=true then animate(implicitplot, [[osc_circle, curve_eq, norm_line], x=x_min..x_max, y=y_min..y_max, numpoints=numPts], t_0=t_begin..t_end, frames=nFrames, scaling=constrained); else animate(implicitplot, [[osc_circle, curve_eq], x=x_min..x_max, y=y_min..y_max, numpoints=numPts], t_0=t_begin..t_end, frames=nFrames, scaling=constrained); end if end proc: |