How can I use procedure, circle3D, to make an animation of the osculating circle of a parametric curve?

Question

I try to use circle3D procedure, shown below, to make an animation of the osculating circle of a parametric curve.

(* circle3D PROCEDURE *)

(* Let's create circle3D that is something we would expect from 
Circle but with an extra argument for its normal vector *)

circle3D[centre_: {0, 0, 0}, radius_: 1, normal_: {0, 0, 1}, angle_: {0, 2 Pi}] :=
  Composition[
    Line,
    Map[RotationTransform[{{0, 0, 1}, normal}, centre], #] &, 
    Map[Append[#, Last@centre] &, #] &,
    Append[DeleteDuplicates[Most@#], Last@#] &,
    Level[#, {-2}] &,
    MeshPrimitives[#, 1] &,
    DiscretizeRegion,
    If][First @ Differences @ angle >= 2 Pi, 
        Circle[Most @ centre, radius],
        Circle[Most @ centre, radius, angle]]

I am a new user of this site and also of Mathematica. I found this procedure on the site. It works perfectly with Graphics3D, but when I try to use it o make an animation I receive an error with the following content:

Coordinate TransformationFunction[{{0.4413884991460323, -0.5698965816793001, 0.693104666707606, -54.82771835995227}, {-0.5698965816793001, 0.4185903560645521, 0.7071067811865474, -55.93534903389071}, {-0.6931046667076058, -0.7071067811865474, -0.14002114478941538`, should be a triple of numbers, or a Scaled form.

Below is the code from my Mathematica notebook:

(* Solution: *)

ClearAll["Global`*"];

(* ZAŁOŻENIE GLOBALNE t > 0 *)
$Assumptions = t > 0;

(* KRZYWA *)
r[t_] := {t + 1/t, t - 1/t, 2 Log[t]};

(* JEDNOSTKOWY WEKTOR STYCZNY *)
tv = D[r[t], t];
tvnorm = Simplify[Sqrt[tv.tv]];
utv = tv/tvnorm;

(* POCHODNA WEKTORA STYCZNEGO *)
tvD = D[tv, t];

(* JEDNOSTKOWY WEKTOR BINORMALNY  *)
bv = tv\[Cross]tvD;
bvnorm = Simplify[Sqrt[bv.bv]];
ubv = bv/bvnorm;

(* JEDNOSTKOWY WEKTOR NORMALNY  *)
unv = ubv\[Cross]utv;

(* ORTOGONALNOŚĆ WEKTORÓW utv, ubv, unv *)
Simplify[(ubv\[Cross]utv)\[Cross]unv];

(* KRZYWIZNA KRZYWEJ *)
k = bvnorm/tvnorm^3;

(* PROMIEŃ OKRĘGU OSKULACYJNEGO *)
or = 1/k;

(* WEKTOR WODZĄCY ŚRODKA OKRĘGU OSKULACYJNEGO *)
occv = r[t] + or*unv;

(* W PUNKCIE M: *)

(* PUNKT M - tutaj t = 1 *)
pM = {2, 0, 0};

(* JEDNOSTKOWY WEKTOR STYCZNY w M *)
utvM = utv /. t -> 1;

(* JEDNOSTKOWY WEKTOR BINORMALNY w M *)
ubvM = ubv /. t -> 1;

(* JEDNOSTKOWY WEKTOR NORMALNY w M *)
unvM = unv /. t -> 1;

(* KRZYWIZNA KRZYWEJ w M *)
kM = k /. t -> 1;

(* PROMIEŃ OKRĘGU OSKULACYJNEGO w M *)
orM = 1/kM;

(* WEKTOR WODZĄCY ŚRODKA OKRĘGU OSKULACYJNEGO w M *)
occvM = r[t] + or*unv /. t -> 1;

In[23]:= (* circle3D PROCEDURE *)

(* Let's create circle3D that is something we would expect from 
   Circle but with an extra argument for its normal vector *)

circle3D[centre_: {0, 0, 0}, radius_: 1, normal_: {0, 0, 1}, 
  angle_: {0, 2 Pi}] :=
 Composition[
   Line,
   Map[RotationTransform[{{0, 0, 1}, normal}, centre], #] &, 
   Map[Append[#, Last@centre] &, #] &,
   Append[DeleteDuplicates[Most@#], Last@#] &,
   Level[#, {-2}] &,
   MeshPrimitives[#, 1] &,
   DiscretizeRegion,
   If
   ][First@Differences@angle >= 2 Pi,
  Circle[Most@centre, radius],
  Circle[Most@centre, radius, angle]
  ]

(* Visualization: *)

In[24]:= (* SKALOWANIE WEKTORÓW JEDNOSTKOWYCH *)
scalvs = 3;

(* KRZYWA *)
curV = ParametricPlot3D[r[t], {t, 0.1, 4}, PlotStyle -> Red];

(* PUNKT M NA KRZYWEJ *)
pMV = Graphics3D[{Blue, PointSize[Large], Point[pM]}];

(* ŚRODEK OKRĘGU OSKULACYJNEGO W M *)
occvMV = Graphics3D[{Black, PointSize[Large], Point[occvM]}];

(* OKRĄG OSKULACYJNY W M *)
ocMV = Graphics3D[circle3D[occvM, orM, ubvM ]];

(* JEDNOSTKOWY WEKTOR STYCZNY W M *)
utvMV = Graphics3D[{Green, Arrow[{pM, pM + scalvs*utvM}]}];

(* JEDNOSTKOWY WEKTOR NORMALNY W M *)
unvMV = Graphics3D[{Magenta, Arrow[{pM, pM + scalvs*unvM}]}];

(* JEDNOSTKOWY WEKTOR BINORMALNY W M *)
ubvMV = Graphics3D[{Cyan, Arrow[{pM, pM + scalvs*ubvM}]}];

(* WYKRES *)
Show[{curV, pMV, occvMV, ocMV, utvMV, ubvMV, unvMV}, 
 BaseStyle -> {FontSize -> 12, FontFamily -> "Verdena"}, Axes -> True, 
 AxesLabel -> {x, y, z}, Ticks -> Automatic, AxesStyle -> {Red, Green, Blue}, 
 PlotRange -> {{1, 12}, {-5, 4}, {-5, 3}}, 
 PlotLabel -> 
  Style[Framed["Krzywa i okrąg oskulacyjny", FrameStyle -> Red], Bold, 14, 
   Black, Background -> Lighter[LightYellow]], 
 BaseStyle -> {FontSize -> 12, FontFamily -> "Verdena"}]        

(* Animation: *)

(* SKALOWANIE WEKTORÓW JEDNOSTKOWYCH *)
scalvs = 3;

(* KRZYWA *)
curV = ParametricPlot3D[r[t], {t, 0.1, 4}, PlotStyle -> Red];

(* PEWIEN PUNKT P NA KRZYWEJ *)
pPV[t_] = Graphics3D[{Blue, PointSize[Large], Point[r[t]]}];

(* ŚRODEK OKRĘGU OSKULACYJNEGO W P *)
occvPV[t_] = Graphics3D[{Black, PointSize[Large], Point[occv]}];

(* OKRĄG OSKULACYJNY W P *)
ocPV[t_] = Graphics3D[circle3D[occv, or, ubv ]];

(* JEDNOSTKOWY WEKTOR STYCZNY W P *)
utvPV[t_] = Graphics3D[{Green, Arrow[{r[t], r[t] + scalvs*utv}]}];

(* JEDNOSTKOWY WEKTOR NORMALNY W P *)
unvPV[t_] = Graphics3D[{Magenta, Arrow[{r[t], r[t] + scalvs*unv}]}];

(* JEDNOSTKOWY WEKTOR BINORMALNY W P *)
ubvPV[t_] = Graphics3D[{Cyan, Arrow[{r[t], r[t] + scalvs*ubv}]}];

(* ANIMATION *)

Animate[
  Show[{curV, pPV[t], occvPV[t], ocPV[t], utvPV[t], ubvPV[t], nvPV[t]},
    BaseStyle -> {FontSize -> 12, 
    FontFamily -> "Verdena"}, 
    Axes -> True, 
    AxesLabel -> {x, y, z}, 
    Ticks -> Automatic, 
    AxesStyle -> {Red, Green, Blue},
    PlotRange -> {{1, 12}, {-5, 4}, {-5, 3}}, 
    PlotLabel -> 
      Style[Framed["Krzywa i okrąg oskulacyjny", 
        FrameStyle -> Red], Bold, 14, Black, 
        Background -> Lighter[LightYellow]], 
    BaseStyle -> {FontSize -> 12, FontFamily -> "Verdena"}], {t, 0.1, 4}, 
    AnimationRunning -> False]

Please help me solve the problem I have described.

Here's an alternative (more efficient) 3D arc procedure that you might find useful: https://mathematica.stackexchange.com/questions/10957/an-efficient-circular-arc-primitive-for-graphics3d — b3m2a1, Dec 02 '17 at 23:15

score 5 · Answer 1 · answered Mar 04 '18 at 02:43

There's no need for the full machinery of Sjoerd's fine answer for this application, since you don't need an arc, but the whole circle. Thus, the simple NURBS representation of a circle, in combination with the use of FrenetSerretSystem[], ought to suffice:

osculatingCircle[fun_?VectorQ, {t_, tvalue_}] := Module[{bv, ka, nv, ta, tv},

          {{ka, ta}, {tv, nv, bv}} = FrenetSerretSystem[fun, t] /. t -> tvalue;

          BSplineCurve[Composition[TranslationTransform[(fun /. t -> tvalue) + nv/ka], 
                                   RotationTransform[{{0, 0, 1}, bv}]][
                       {{1, 0, 0}, {1, 1, 0}, {-1, 1, 0}, {-1, 0, 0},
                        {-1, -1, 0}, {1, -1, 0}, {1, 0, 0}}/ka], 
                       SplineDegree -> 2, 
                       SplineKnots -> {0, 0, 0, 1/4, 1/2, 1/2, 3/4, 1, 1, 1}, 
                       SplineWeights -> {1, 1/2, 1/2, 1, 1/2, 1/2, 1}]]

As an example:

tref = KnotData["Trefoil", "SpaceCurve"];

Animate[Show[ParametricPlot3D[tref[t], {t, 0, 2 π}], 
             Graphics3D[{Directive[Red, AbsoluteThickness[2]], 
                         osculatingCircle[tref[t], {t, N[tm, 20]}]}], 
             PlotRange -> 10], {tm, 0, 2 π, π/20}]

trefoil and its osculating circle

So compact! I never knew FrenetSerretSystem existed (I never would have even guessed that was called the Frenet-Serret basis). That seems very useful. — b3m2a1, Mar 04 '18 at 03:29
It is very handy. I remember having to code all that functionality manually 15 years ago, so I am glad it's built-in now. — J. M.'s missing motivation, Mar 04 '18 at 12:06

b3m2a1 · Answer 2 · 2018-01-02T16:11:12.957

Your issue is the circle3D procedure as you can see from trying this:

With[{t = 4},
 Show[{curV, pPV[t], occvPV[t],(*ocPV[t],*)utvPV[t], ubvPV[t], 
   unvPV[t]}, BaseStyle -> {FontSize -> 12, FontFamily -> "Verdena"}, 
  Axes -> True, AxesLabel -> {x, y, z}, Ticks -> Automatic, 
  AxesStyle -> {Red, Green, Blue}, 
  PlotRange -> {{1, 12}, {-5, 4}, {-5, 3}}, 
  PlotLabel -> 
   Style[Framed["Krzywa i okrąg oskulacyjny", FrameStyle -> Red], 
    Bold, 14, Black, Background -> Lighter[LightYellow]], 
  BaseStyle -> {FontSize -> 12, FontFamily -> "Verdena"}]
 ]

I'd suggest adapting the method here to work with your circle3D input.

Perhaps like so:

ClearAll[sjoerdCircle];
sjoerdCircle[m_List, r_, angles_List: {0, 2 \[Pi]}] := 
  Module[{seg, \[Phi], start, end, pts, w, k}, {start, end} = 
     Mod[angles // N, 2 \[Pi]];
    If[end <= start, end += 2 \[Pi]];
    seg = Quotient[end - start // N, \[Pi]/2];
    \[Phi] = Mod[end - start // N, \[Pi]/2];
    If[seg == 4, seg = 3; \[Phi] = \[Pi]/2];
    pts = 
     r RotationMatrix[start].# & /@ 
      Join[Take[{{1, 0}, {1, 1}, {0, 1}, {-1, 1}, {-1, 
          0}, {-1, -1}, {0, -1}}, 2 seg + 1], 
       RotationMatrix[seg \[Pi]/2].# & /@ {{1, 
          Tan[\[Phi]/2]}, {Cos[\[Phi]], Sin[\[Phi]]}}];
    If[Length[m] == 2, pts = m + # & /@ pts, 
     pts = m + # & /@ 
       Transpose[
        Append[Transpose[pts], ConstantArray[0, Length[pts]]]]];
    w = Join[
      Take[{1, 1/Sqrt[2], 1, 1/Sqrt[2], 1, 1/Sqrt[2], 1}, 
       2 seg + 1], {Cos[\[Phi]/2], 1}];
    k = Join[{0, 0, 0}, Riffle[#, #] &@Range[seg + 1], {seg + 1}];
    BSplineCurve[pts, SplineDegree -> 2, SplineKnots -> k, 
     SplineWeights -> w]] /; Length[m] == 2 || Length[m] == 3;
sjoerdCircleNormal[
  centre : {Repeated[_?NumericQ, {3}]} : {0, 0, 0},
  radius : _?NumericQ : 1,
  normal : {Repeated[_?NumericQ, {3}]} : {0, 0, 1},
  angle : {Repeated[_?NumericQ, {2}]} : {0, 2 Pi}
  ] :=
 With[{norm = Normalize[normal]},
  Which[
   norm == {0, 0, 1},
   sjoerdCircle[centre, radius, angle],
   norm == {0, 0, -1},
   sjoerdCircle[centre, radius, \[Pi] + angle],
   True,
   GeometricTransformation[
    sjoerdCircle[centre, radius, angle],
    RotationTransform[{{0, 0, 1}, normal}, centre]
    ]
   ]
  ]

Then you can do this:

ocPV[t_] = Graphics3D[sjoerdCircleNormal[occv, or, ubv]];

And your animation wil work fine.

Thank you for your answer. Animation works. However, the call ocPV[t_] = Graphics3D [sjoerdCircleNormal [occv, or, ubv]]; in the animation does not create an osculation circle, but rather a circle that intersects the curve, and the center of such an animated circle does not coincide with occv. I do not know how to improve it. Can you help me? — SIJA, Dec 03 '17 at 11:25
@SIJA you just needed to recenter the RotationTransform. It does what you want now. — b3m2a1, Jan 02 '18 at 16:12

score 1 · Answer 3 · answered Nov 26 '20 at 14:51

My Manipulate version:

OsculatingCircle[r_, tr_, opt : OptionsPattern[]] := 
 DynamicModule[{\[Kappa], Tv, Nv, Cv, rc, s1, t1, t2},
  t1 = tr[[1]];
  \[Kappa] = 
   ComplexExpand[Norm[D[r, t1]\[Cross]D[r, {t1, 2}]]/Norm[D[r, t1]]^3];
  Tv = ComplexExpand[D[r, t1]/Norm[D[r, t1]]];
  Nv = ComplexExpand[D[Tv, t1]/Norm[D[Tv, t1]]];
  Cv = Nv/\[Kappa] + r;
  rc = Cv + Sin[s1] Nv/\[Kappa] + Cos[s1] Tv/\[Kappa];
  Manipulate[
   Show[ParametricPlot3D[r, tr], 
    ParametricPlot3D[Evaluate[rc /. tr[[1]] -> t2], {s1, 0, 2 Pi}, 
     PlotStyle -> Orange]], {{t2, tr[[2]], tr[[1]]}, tr[[2]], 
    tr[[3]]}, ControlPlacement -> Bottom]
  ]

For example:

OsculatingCircle[{t - 3/2 Sin[t], 1 - 3/2 Cos[t], t}, {t, 0, 6 Pi}]

How can I use procedure, circle3D, to make an animation of the osculating circle of a parametric curve?

3 Answers3

Linked