import numpy as np
from scipy.special import ellipj, ellipk
import ddg
[docs]def jacobi_elliptic_curve(k):
"""Returns the elliptic curve [cn, sn, dn] for a parameter k.
More precisely, returns the collection of curves obtained by first taking
the curve ::
u -> [cn(u, k), sn(u, k), dn(u, k)]
and adding the reflection of this about the plane z = 0. In the case where
k = 1, we also have to reflect about the planes x = 0 and x + z = 0.
This curve (collection) is the intersection of the two cylinders ::
x^2 + y^2 = 1,
k^2 * y^2 + z^2 = 1
Homogenizing in the last component, to [cn, sn, dn, 1], gives the
intersection curve of any two quadrics in the pencil ::
u1 * diag([1, 1, 0, -1]) + u2 * diag([0, k, 1, -1])
in homogeneous coordinates.
Parameters
----------
k : float
Parameter in the interval [0, 1]. It is best to avoid values close
(about 1e-9) but not equal to 1 because the scipy.ellipj functions do
not work well for those values.
Returns
-------
NetCollection
"""
def curve_fct(u):
sn, cn, dn, ph = ellipj(u, k)
return np.array([cn, sn, dn])
# inf if k == 1
K = ellipk(k)
if K < np.inf:
domain = [[-2 * K, 2 * K, True]]
else:
domain = [[-K, K]]
curves = [ddg.SmoothNet(curve_fct, domain)]
trafos = [np.diag([1, 1, -1])]
if k == 1:
trafos += [np.diag([-1, 1, 1]), np.diag([-1, 1, -1])]
for trafo in trafos:
curve = ddg.SmoothNet(curve_fct, domain)
curve.transform(trafo)
curves.append(curve)
return ddg.NetCollection(curves)