This animation uses the modular flow to produce a periodic family of Jacobi elliptic functions of type cn, visualized using domain coloring. This is also explained in my paper "Looping animations using the modular flow and elliptic functions."
ABSTRACT (original article): This paper describes an approach to generating looping animations using the modular flow and elliptic functions. The modular flow is a flow on lattices with many periodic orbits, and elliptic functions are meromorphic, doubly-periodic functions which can be visualized using domain coloring. CITATION (original article): Clayton Shonkwiler (2026), Looping Animations Using the Modular Flow and Elliptic Functions, Proceedings of Bridges 2026: Mathematics and the Arts, 533--536. http://archive.bridgesmathart.org/2026/bridges2026-533.html

The Modular Flow

See Quasi or the paper for an explanation of the modular flow. The upshot is that when
Λ
is the planar lattice generated by
u
1
=
ϕ
1

and
u
2
=
1
-ϕ

(
ϕ=
1+
5
2
is the golden ratio), and we act on the lattice points by the family of matrices
t

0
0
-t

to get the family
t·Λ
of lattices, this family turns out to be periodic with period
t
0
=ln(1+ϕ)
:
In[]:=
DynamicModule[{M,pts,U={{GoldenRatio,1},{-1,GoldenRatio}}},​​pts=Flatten[Table[{x,y},{x,-15,15},{y,-15,15}],1];​​Animate[​​M=DiagonalMatrix[{
s
E
,
-s
E
}].U;​​Graphics[{Thickness[.005],PointSize[Large],Point[M.#]&/@pts},PlotRange->10],​​{s,0,Log[1+GoldenRatio]},SaveDefinitions->True,AnimationRunning->False]​​]
Out[]=
s
In particular, if
U
has columns
u
1
and
u
2
and
B=
2
1
1
1

, then
diag(
t
0

,
-
t
0

)U=UB
or, equivalently,
diag(
t
0

,
-
t
0

)
u
1
=2
u
1
+
u
2
and
diag(
t
0

,
-
t
0

)
u
2
=
u
1
+
u
2
.

The Jacobi elliptic function

In general, elliptic functions are doubly-periodic meromorphic functions on the complex plane

. The Jacobi elliptic functions are a basic family of elliptic functions of the form
pq(u|m)
, where p and q are any of the letters c, s, n, or d,
u
is the argument, and
m
is a parameter which implicitly specifies the lattice (or, equivalently, a given lattice specifies the parameter
m
): for any integers
α
and
β
, ​
pq(u+2αK(m)+2βK(1-m)|m)=
γ
(-1)
pq(u|m)
, ​where
γ
is determined by p, q,
α
, and
β
and
K
is the elliptic integral of the first kind. For example, for the cn function we have
γ=α+β
, so that ​
cn(u+4αK(m)|m)=cn(u|m)
and
cn(u+4βK(1-m)|m)=cn(u|m)
, and we see that the periods of cn are
4K(m)
and
4K(1-m)
.​For example,
cn
is implemented in Mathematica as JacobiCN:
In[]:=
ComplexPlot[JacobiCN[u,.8],{u,10}]
Out[]=
In general,
pq(u|m)
will have a zero in the p corner and a pole in the q corner of the rectangle (or, more generally, parallelogram) coming from the lattice, where the corners are labeled as in this Wikipedia diagram: ​​
​​For this parameter
m=0.8
, we can compute
K(m)
and
K(1-m)
as follows:
In[]:=
EllipticK[0.8]
Out[]=
2.25721
In[]:=
EllipticK[1-0.8]
Out[]=
1.65962
So we expect to see a (simple) zero of
cn(·|0.8)
at
≈2.257
and a (simple) pole at
≈1.659
, which is indeed what we see in the plot above.
Now, the idea is that we’re going to get an animation by plotting the family
cn(·|
m
t
)
of Jacobi functions, where
m
t
is the family of parameters corresponding to our family
t·Λ
of lattices. In principle, since we know the periods of cn and we know the lattice
t·Λ
, we could try to solve for
m
t
. I don’t doubt that some cleverer person than I could figure this out, but I spent a very long time trying and failing to figure this out due to the facts that: (i)
K
has a branch cut discontinuity running from 1 to ∞ and (ii) I’m representing
t·Λ
by
t·
u
1
=
t

ϕ
-t

and
t·
u
2
=
t

-
-t

ϕ
, which span an increasingly long and pointy parallelogram.
In[]:=
AnimateGraphicsParallelogram[{0,0},{{
t
E
GoldenRatio,
-t
E
},{
t
E
,-
-t
E
GoldenRatio}}],PlotRange->5,{t,0,Log[1+GoldenRatio]},SaveDefinitions->True,AnimationRunning->False
Out[]=
t
So instead I’m going to represent cn in terms of the Weierstrass
℘
-function. In general, an even elliptic function
f
with zeros at
±
z
1
,…,±
z
n
and poles at
±
p
1
,…,±
p
n
can be represented as​
f(z)=a
n
∏
i=1
℘(z)-℘(
z
i
)
℘(z)-℘(
p
i
)
, where
a
is some fixed constant (see, for example, §8.10 in Lawden’s Elliptic Functions and Applications). Since cn is even, we can use this representation. We know from the construction of cn that it should have zeros at
±
ω
1
4
and
±
3
ω
1
4
+
ω
2
2
and poles at
±
ω
2
4
and
±
3
ω
2
4
+
ω
1
2
, so we can just plug these into the above representation to get an alternative formula for cn.
Translating to complex numbers, we saw above that after one period
In particular, our zeros and poles get sent as follows:
So we’ll need to run the modular flow for 3 periods to get a looping animation.
The last piece is to use the vik0 cyclic gradient from Scientific Color Maps:
So here’s the code for the animation:
Which we can export to a GIF:
Or visualize using ListAnimate:
To get the final GIF under 10 MB, I used Gifsicle, specifically the following command:

CITE THIS NOTEBOOK

cn (elliptic function animation)​
by Clayton Shonkwiler​
Wolfram Community, STAFF PICKS, June 24, 2026
​https://community.wolfram.com/groups/-/m/t/3738945