1
$\begingroup$

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.

Improve this question
$\endgroup$
1

3 Answers 3

Reset to default
5
$\begingroup$

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

Improve this answer
$\endgroup$
2
  • 1
    $\begingroup$ So compact! I never knew FrenetSerretSystem existed (I never would have even guessed that was called the Frenet-Serret basis). That seems very useful. $\endgroup$
    – b3m2a1
    Commented Mar 4, 2018 at 3:29
  • $\begingroup$ 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. $\endgroup$
    – J. M.'s missing motivation ♦
    Commented Mar 4, 2018 at 12:06
1
$\begingroup$

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.

Improve this answer
$\endgroup$
2
  • $\begingroup$ 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? $\endgroup$
    – SIJA
    Commented Dec 3, 2017 at 11:25
  • $\begingroup$ @SIJA you just needed to recenter the RotationTransform. It does what you want now. $\endgroup$
    – b3m2a1
    Commented Jan 2, 2018 at 16:12
1
$\begingroup$

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}]

enter image description here

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.