# Make a "formula-millipede" move along arbitrary path

Math artist Chirag Dudhat specializes in simple math that does something interesting visually. He recently had a hit with a "millipede" simulation with simple complex plane formulas. Below is a small part and full video is here.

The full video has a frame with formulas he used:

Quoting Chirag Dudhat own words:

The equation displayed in the video , represents the movement of a single leg. Zroot is where the leg connects to the body, Zjoint is the point where the leg bends and Ztip is the free point of the leg which may touch the surface. Change the value of L from 1 to n for n number of legs. For example , L = 12 gives you the position of 12th leg on complex plane and changing the value of "t" will animate it. This equation gives you a straight walking millipede. However for this animation I drew a curved path on my iPad , exported svg file, extracted x and y coordinates of the path. Moving the centipede on this path is like attaching the legs to the path itself and sliding the whole structure on the curve.

But it turned out to be a bit tricky to get from simple to final form. Here is what I tried. First of all his formulas have a bit different parameter numbers to get the effect and furthermore - another formula is needed to get decreasing size of legs, thier complete number, and some other details. So here is my set of formulas:

maxL=73;
zR[L_]:=0.5I L

damp[L_,max_]:=Exp[-((L-max/2)^2/(15max))^3]

zJ[L_,t_]:=zR[L]+(-1)^L 3Exp[(5π I)/180 Sin[t+(π L)/6]]damp[L,maxL]
zT[L_,t_]:=zR[L]+(-1)^L 5Exp[(17π I)/180 Sin[t+(π L)/6]]damp[L,maxL]

zJs[L_,t_]:=zR[L]+(-1)^(L+1) 3Exp[(5π I)/180 Sin[t+(π L)/6]]damp[L,maxL]
zTs[L_,t_]:=zR[L]+(-1)^(L+1) 5Exp[(17π I)/180 Sin[t+(π L)/6]]damp[L,maxL]


This is easily animated with a Manipulate:

Manipulate[Rotate[Graphics[{
Point/@ReIm@Table[{zR[L],zJ[L,t],zT[L,t],zJs[L,t],zTs[L,t]},{L,maxL}],
Line[ReIm@Table[zR[L],{L,maxL}]],
Table[Line/@ReIm@{{zR[L],zJ[L,t],zT[L,t]},{zR[L],zJs[L,t],zTs[L,t]}},{L,maxL}]
},ImageSize->{200,600}],Pi/2],{t,0,33},SaveDefinitions->True]


## My question:

How can we make this millipede to move along a path in the plane? For instance a closed random loop given by spline or any other similar path:

spline=BSplineFunction[RandomReal[1,{10,2}],SplineClosed->True]
ParametricPlot[spline[t],{t,0,1}]


Also, can we simplify all the code somehow? Feel free to get creative and create your own version if you want to. All millipedes moving along a curved path are cool :-)

• The segments of the "millipede" is segments=60.
• We using ArcLength since the length of millipede always a constant.
• The profile of the "millipede" seems like a ellipsoid, so we try to use the ellise curve equation of λ = 3 Sqrt[1 - (L - segments/2)^2/(segments/2)^2].
• e1,e2 are the moving orthonormal frame.
• Maybe AnglePath can simplify the code partly,but we need times to figure out.
• We ajust the angular velocity  ω to .2.
Clear["Global*"];
pts = {{0, 0}, {1, 1}, {2, -1}, {3, 0}, {4, -2}, {6, 1}, {2, 3}, {2,
0}, {3, -2}, {4, -1}, {5, 2}, {4, 2}, {3, 2}, {1, 1}, {1, -1}};

f = BSplineFunction[pts];
segments = 60;
factor1 = 4;
factor2 = 3;
n = Floor[factor1*segments];
plot0 = ParametricPlot[f[t], {t, 0, 1},
MeshFunctions -> {"ArcLength"}, Mesh -> n, MaxRecursion -> 0,
PlotPoints -> Floor[factor1*factor2*segments]];
positions = Cases[Normal@plot0, Point[pts_] :> pts, -1];
frame[t_, L_] :=
Module[{center, nest, vector, norm, e2, e1, λ, ω,
dirs},
center = positions[[t + L]];
nest = positions[[t + L + 1]];
vector = nest - center;
norm = Norm[vector];
e2 = Normalize[vector];
e1 = -Cross@e2;
λ = 3*Sqrt[1 - (L - segments/2)^2/(segments/2)^2];
ω = .2;
dirs = {AngleVector[{norm*λ,
2*5*Degree*Sin[ω*t + π/6*L]}],
AngleVector[{norm*λ,
2*(5 + 17)*Degree*Sin[ω*t + π/6*L]}]} . {e1, e2};
Graphics[{White, Line[{center, nest}], AbsolutePointSize[2],
AbsoluteThickness[1],
Through@{Point, Line}@FoldList[Plus, center, -dirs],
Through@{Point, Line}@FoldList[Plus, center, dirs]}]];
list = Table[
Show[Table[frame[t, L], {L, 1, segments - 1}],
AspectRatio -> Automatic, Background -> Black, PlotRange -> All,
ImageSize -> {400, 400}], {t, 1, n - segments - 1, 1}];
ani = ListAnimate[list, AnimationRate -> 10]


• Test pts = Table[{Sin@s, s}, {s, 0, 8 π, .1}];

• We can use Canvas[] to get some arbitrary path and smoothing it by MovingAverage.

img;
pts = img[[1, 1]];
pts = pts[[1 ;; -1 ;; 10]];
pts = MovingAverage[pts, 20];


## Try to simulate the lighting

by using RadialGradientFilling[ColorData["Gradients"][[11]]] or RadialGradientFilling[ColorData["Gradients"][[24]]]

Clear["Global*"];
pts = {{0, 0}, {1, 1}, {2, -1}, {3, 0}, {4, -2}, {6, 1}, {2, 3}, {2,
0}, {3, -2}, {4, -1}, {5, 2}, {4, 2}, {3, 2}, {1, 1}, {1, -1}};

f = BSplineFunction[pts];
segments = 60;
factor1 = 5;
factor2 = 3;
n = Floor[factor1*segments];
plot0 = ParametricPlot[f[t], {t, 0, 1},
MeshFunctions -> {"ArcLength"}, Mesh -> n, MaxRecursion -> 0,
PlotPoints -> Floor[factor1*factor2*segments]];
positions = Cases[Normal@plot0, Point[pts_] :> pts, -1];
frame[t_, L_] :=
Module[{center, nest, vector, norm, e2, e1, λ, ω,
dirs, reg},
center = positions[[t + L]];
nest = positions[[t + L + 1]];
vector = nest - center;
norm = Norm[vector];
e2 = Normalize[vector];
e1 = -Cross@e2;
λ = 3*Sqrt[1 - (L - segments/2)^2/(segments/2)^2];
ω = .2;
dirs = {AngleVector[{norm*λ,
2*5*Degree*Sin[ω*t + π/6*L]}],
AngleVector[{norm*λ,
2*(5 + 17)*Degree*Sin[ω*t + π/6*L]}]} . {e1, e2};
Partition[FoldList[Plus, center, -dirs], 2, 1],
Partition[FoldList[Plus, center, dirs], 2, 1],
Disk[#, norm/2] & /@ FoldList[Plus, center, -dirs],
Disk[#, norm/2] & /@ FoldList[Plus, center, dirs]}];
ani = Animate[
Table[frame[t, L], {L, 1, segments - 1}]}, Background -> Black,
ImageSize -> {400, 400}, PlotRange -> All], {t, 1, n - segments,
1}, AnimationRate -> 10]


• +1 thank you. I see you are trying to post a nice GiF. Something that helped me... One can go HERE and use THIS in the sequence shown. It will highly compress your GIF and allow for longer / smoother frames if you want to. Commented Jan 8 at 6:41
• @VitaliyKaurov Thanks! Commented Jan 8 at 10:37

Change zR to add a wiggle, choose a path g, and add a transform function f:

Clear[maxL, zR, damp, zJ, zT, zJs, zTs, f, g];

maxL = 73;
zR[L_, t_] :=(**)-Exp[-2 (5 \[Pi] I)/
180 Sin[-t + (\[Pi] L)/24]](**).5 I L

g[t_] := {Sin[(Pi t/33)], Cos[(Pi t/33)]};
f[q_, t_] :=
11 g[t] +
RotationMatrix[
Quiet[ArcTan[Evaluate@D[g[x], x][[1]],
Evaluate@D[g[x], x][[2]]] /. {Indeterminate -> 0, x -> t}]] .
RotationMatrix[-Pi/2] . q;

damp[L_, max_] := Exp[-((L - max/2)^2/(15 max))^3]

zJ[L_, t_] :=
zR[L, t] + (-1)^L 3 Exp[(5 \[Pi] I)/180 Sin[t + (\[Pi] L)/6]] damp[L,
maxL]
zT[L_, t_] :=
zR[L, t] + (-1)^L 5 Exp[(17 \[Pi] I)/180 Sin[t + (\[Pi] L)/6]] damp[
L, maxL]

zJs[L_, t_] :=
zR[L, t] + (-1)^(L + 1) 3 Exp[(5 \[Pi] I)/180 Sin[
t + (\[Pi] L)/6]] damp[L, maxL]
zTs[L_, t_] :=
zR[L, t] + (-1)^(L + 1) 5 Exp[(17 \[Pi] I)/180 Sin[
t + (\[Pi] L)/6]] damp[L, maxL]

Manipulate[
Rotate[Graphics[{Point /@ (f[#, t] & /@ # & /@
ReIm@Table[{zR[L, t], zJ[L, t], zT[L, t], zJs[L, t],
zTs[L, t]}, {L, maxL}]),
Line /@ Table[
f[#, t] & /@ # & /@
ReIm@{{zR[L, t], zJ[L, t], zT[L, t]}, {zR[L, t], zJs[L, t],
zTs[L, t]}}, {L, maxL}]}, ImageSize -> {600, 600},
PlotRange -> {{-#, #}, {-#, #}} &[40]], Pi/2], {t, 1, 66,
AnimationRate -> 1, AnimationRunning -> False,
Appearance -> "Open"}, SaveDefinitions -> True]


A different approach:

Clear[g, f, zR, damp, zJ, zT, zJs, zTs];
g[t_] :=(*1.5{Sin[3t]+Sin[6(3Sin[3t]/2)]/30,Sin[2 t]+Sin[6(3Sin[2 t]/
2)]/30};*){3 Sin[3 t]/2, 3 Sin[2 t]/2};
f[q_, t_] :=
g[t] + RotationMatrix[
Quiet[ArcTan[Evaluate@D[g[x], x][[1]],
Evaluate@D[g[x], x][[2]]] /. {Indeterminate -> 0, x -> t}]] .
RotationMatrix[-Pi/2] . q;
zR[L_, t_] := 0.5 I L; maxL = 73;
fr = maxL; wormfatness = 1.5; c1 = 6; c2 = 10; c3 = 10; c4 = 17; c5 = \
12;
damp[L_, max_] := Exp[-((L - max/2)^2/(15 max))^3];
zJ[L_, t_] :=
zR[L, t] - (-1)^
L c1 Exp[(c3 \[Pi] I)/180 Sin[fr t + (\[Pi] L)/c5]] damp[L, maxL];
zT[L_, t_] :=
zR[L, t] - (-1)^
L c2 Exp[(c4 \[Pi] I)/180 Sin[fr t + (\[Pi] L)/c5]] damp[L, maxL];
zJs[L_, t_] :=
zR[L, t] - (-1)^(L + 1) c1 Exp[(c3 \[Pi] I)/180 Sin[
fr t + (\[Pi] L)/c5]] damp[L, maxL];
zTs[L_, t_] :=
zR[L, t] - (-1)^(L + 1) c2 Exp[(c4 \[Pi] I)/180 Sin[
fr t + (\[Pi] L)/c5]] damp[L, maxL];
q[func_, t_, L_] := {wormfatness #, #2 - maxL/4} & @@
ReIm@func[L maxL + maxL/2, t];
Manipulate[
Show[(*ParametricPlot[g[x],{x,0 ,2Pi}],*)
Graphics[{White, Thickness[.002], PointSize[.007],
Line[Table[
f[#/maxL, t + L/2] &@q[zR, t, L], {L, -.5, .5, 1/maxL}]],
Table[
Point[f[#/maxL, t + L/2] & /@ {q[zJ, t, L], q[zT, t, L],
q[zJs, t, L], q[zTs, t, L]}], {L, -.5, .5, 1/maxL}],
Table[
Line[f[#/maxL, t + L/2] & /@ {q[zT, t, L],
q[zTs, t, L]}], {L, -.5, .5, 1/maxL}]
}, Background -> Black,
PlotRange -> {{-#, #}, {-#, #}} &[2]]], {t, 0, 2 Pi,
AnimationRate -> 1/(15 Pi), AnimationRunning -> False,
Appearance -> "Open"}, SaveDefinitions -> True]


• +1 thank you. Just FYI about GiFs in case your need it. Something that helped me... One can go here HERE and use THIS in the sequence shown. It will highly compress your GIF and allow for longer / smoother frames if you want to. YOu might not need it though - your gifs look great :-) Commented Jan 8 at 6:42
• Brilliant, thank you, that is useful! Commented Jan 8 at 7:16