"""Coordinate functions for smooth and discrete confocal conics and quadrics."""
import numpy as np
from scipy.special import gamma
from ddg.math.functions import cn, dn, qgamma, sn
sqrt = np.math.sqrt
cos = np.math.cos
sin = np.math.sin
tan = np.math.tan
cosh = np.math.cosh
sinh = np.math.sinh
tanh = np.math.tanh
exp = np.math.exp
# Jacobi elliptic functions
# import sympy
# def sn(u, k):
# return sympy.mpmath.ellipfun('sn', u, k=k)
# def cn(u, k):
# return sympy.mpmath.ellipfun('cn', u, k=k)
# def dn(u, k):
# return sympy.mpmath.ellipfun('dn', u, k=k)
[docs]def confocal_conics_sqrt(u1, u2, a1, a2):
"""confocal conics -- elliptic coordinates (parametrization with square roots)"""
x1_sign = 1
x2_sign = 1
u1 = ((u1 + a1) % (4 * (a1 - a2))) - a1
if -a2 <= u1 < a1 - 2 * a2:
u1 = -u1 - 2 * a2
x2_sign *= -1
elif a1 - 2 * a2 <= u1 < 2 * a1 - 3 * a2:
u1 = u1 - 2 * (a1 - a2)
x1_sign *= -1
x2_sign *= -1
elif 2 * a1 - 3 * a2 <= u1 < 3 * a1 - 4 * a2:
u1 = -u1 + 2 * a1 - 4 * a2
x1_sign *= -1
if u2 < -a2:
u2 = -u2 - 2 * a2
x2_sign *= -1
x1 = sqrt(u1 + a1) * sqrt(u2 + a1) / sqrt(a1 - a2)
x1 *= x1_sign
x2 = sqrt(-(u1 + a2)) * sqrt(u2 + a2) / sqrt(a1 - a2)
x2 *= x2_sign
return np.array([x1, x2])
[docs]def confocal_conics_trigonometric(u1, u2, a1, a2):
"""confocal conics -- parametrization with trigonometric functions"""
def f1(u1, u2, a1, a2):
return sqrt(a1 - a2) * cos(u1) * cosh(u2)
def f2(u1, u2, a1, a2):
return sqrt(a1 - a2) * sin(u1) * sinh(u2)
return np.array([f1(u1, u2, a1, a2), f2(u1, u2, a1, a2)])
[docs]def confocal_conics_concentric(u1, u2, a1, a2):
"""confocal conics -- parametrization with concetric circles on diagonals"""
def f1(u1, u2, a1, a2):
return u1 * u2
def f2(u1, u2, a1, a2):
return sqrt(1 - u1**2) * sqrt(u2**2 - 1)
return np.array([f1(u1, u2, a1, a2), f2(u1, u2, a1, a2)])
[docs]def confocal_conics_ic(u1, u2, a1, a2, k):
"""confocal conics -- parametrization with elliptic coordinates with
diagonals on lines"""
def f1(u1, u2, a1, a2, k):
return sqrt(a1 - a2) * (sn(u1, k) * dn(u2, k)) / (k * cn(u2, k))
def f2(u1, u2, a1, a2, k):
return sqrt(a1 - a2) * (sqrt(1 - k**2) * cn(u1, k)) / (k * cn(u2, k))
return np.array([f1(u1, u2, a1, a2, k), f2(u1, u2, a1, a2, k)])
[docs]def confocal_conics_ic_hyperbolic(u1, u2, a1, a2, k):
"""confocal conics -- parametrization with elliptic coordinates with
diagonals on lines"""
def f1(u1, u2, a1, a2, k):
return k * sqrt(a1 - a2) * sn(u1, k) / sn(u2, k)
def f2(u1, u2, a1, a2, k):
return sqrt(a1 - a2) * dn(u1, k) * cn(u2, k) / sn(u2, k)
return np.array([f1(u1, u2, a1, a2, k), f2(u1, u2, a1, a2, k)])
[docs]def confocal_conics_hyperbolic_pencil(u1, u2):
"""confocal conics -- parametrization st diagonals are vertical lines and circles"""
def f1(u1, u2):
return exp(u1 + u2)
def f2(u1, u2):
return sqrt((1 - exp(2 * u1)) * (exp(2 * u2) - 1))
return np.array([f1(u1, u2), f2(u1, u2)])
[docs]def discrete_confocal_conics_gamma(m, n, a1, a2):
"""discrete confocal conics -- parametrization with gamma functions"""
x1_sign = 1
x2_sign = 1
m = ((m + a1) % (4 * (a1 - a2))) - a1
if -a2 <= m < a1 - 2 * a2:
m = -m - 2 * a2
x2_sign *= -1
elif a1 - 2 * a2 <= m < 2 * a1 - 3 * a2:
m = m - 2 * (a1 - a2)
x1_sign *= -1
x2_sign *= -1
elif 2 * a1 - 3 * a2 <= m < 3 * a1 - 4 * a2:
m = -m + 2 * a1 - 4 * a2
x1_sign *= -1
if n < -a2:
n = -n - 2 * a2
x2_sign *= -1
def D(a1, a2):
return 1 / sqrt(a1 - a2 - 0.5)
def f1(n1, n2, a1, a2):
return (
x1_sign
* D(a1, a2)
* (gamma(n1 + a1 + 0.5) / gamma(n1 + a1))
* (gamma(n2 + a1) / gamma(n2 + a1 - 0.5))
)
def f2(n1, n2, a1, a2):
return (
x2_sign
* D(a1, a2)
* (gamma(-n1 - a2 + 0.5) / gamma(-n1 - a2))
* (gamma(n2 + a2 + 0.5) / gamma(n2 + a2))
)
return np.array([f1(m, n, a1, a2), f2(m, n, a1, a2)])
[docs]def discrete_confocal_conics_gamma_u1_of_corresponding_continuous(n1):
return n1 - 0.5
[docs]def discrete_confocal_conics_gamma_u2_of_corresponding_continuous(n2):
return n2 - 1.0
[docs]def discrete_confocal_conics_trigonometric(m, n, a1, a2, h1, h2, d1=0.0, d2=0.0):
"""discrete confocal conics -- parametrization with trigonometric functions"""
# TODO: Check d1 d2 dependence!!! And add this dependence to the
# 'corresponding_continuous'.
def A(a1, a2, h1, h2):
return sqrt((a1 - a2) / (cos(h1 / 2) * cosh(h2 / 2)))
def x1(n1, n2, a1, a2, h1, h2, d1, d2):
return A(a1, a2, h1, h2) * cos(h1 * n1 + d1) * cosh(h2 * n2 + d2)
def x2(n1, n2, a1, a2, h1, h2, d1, d2):
return A(a1, a2, h1, h2) * sin(h1 * n1 + d1) * sinh(h2 * n2 + d2)
return np.array(
[x1(m, n, a1, a2, h1, h2, d1, d2), x2(m, n, a1, a2, h1, h2, d1, d2)]
)
[docs]def discrete_confocal_conics_trigonometric_u1_of_corresponding_continuous(
n1, a1, a2, h1, d1=0.0
):
"""Returns u1 of the corresponding continuous conics in the sqrt parametrization.
The continuous conics corresponding to (u1, u2) in the sqrt parametrization
establish the polarity relation between points and edges of two dual meshes of
the discrete_confocal_conics_trigonometric parametrization.
"""
return (
-a1 * sin(h1 * n1 + d1) ** 2
- a2 * cos(h1 * n1 + d1) ** 2
- (a1 - a2) * sin(h1 * n1 + d1) * cos(h1 * n1 + d1) * tan(h1 / 2)
)
[docs]def discrete_confocal_conics_trigonometric_u2_of_corresponding_continuous(
n2, a1, a2, h2, d2=0.0
):
"""Returns u1 of the corresponding continuous conics in the sqrt parametrization.
The continuous conics corresponding to (u1, u2) in the sqrt parametrization
establish the polarity relation between points and edges of two dual meshes of
the discrete_confocal_conics_trigonometric parametrization.
"""
return (
a1 * sinh(h2 * n2 + d2) ** 2
- a2 * cosh(h2 * n2 + d2) ** 2
+ (a1 - a2) * sinh(h2 * n2 + d2) * cosh(h2 * n2 + d2) * tanh(h2 / 2)
)
[docs]def discrete_confocal_conics_concentric(m, n, a1, a2, d, c1=0.0, c2=0.0):
"""discrete confocal conics -- diagonally related to concetric circles"""
def x1(n1, n2, a1, a2, d, c1, c2):
return sqrt(a1 - a2) * d**2 * (n1 + c1) * (n2 + c2)
def g1(n1, d, c1):
n1m = -c1 - 1 / 4 - sqrt(16 + d**2) / (4 * d)
n1p = -c1 - 1 / 4 + sqrt(16 + d**2) / (4 * d)
return (
d
* (gamma(n1 - n1m + 1 / 2) * gamma(-n1 + n1p + 1.0))
/ (gamma(n1 - n1m) * gamma(-n1 + n1p + 1 / 2))
)
def g2(n2, d, c2):
n2m = -c2 - 1 / 4 - sqrt(16 + d**2) / (4 * d)
n2p = -c2 - 1 / 4 + sqrt(16 + d**2) / (4 * d)
return (
d
* (gamma(n2 - n2m + 1 / 2) * gamma(n2 - n2p + 1 / 2))
/ (gamma(n2 - n2m) * gamma(n2 - n2p))
)
def x2(n1, n2, a1, a2, d, c1, c2):
return sqrt(a1 - a2) * g1(n1, d, c1) * g2(n2, d, c2)
return np.array([x1(m, n, a1, a2, d, c1, c2), x2(m, n, a1, a2, d, c1, c2)])
[docs]def discrete_confocal_conics_concentric_u1_of_corresponding_continuous(
n1, a1, a2, d, c1, c2
):
return (a1 - a2) * d**2 * (n1 + c1) * (n1 + c1 + 1 / 2) - a1
[docs]def discrete_confocal_conics_concentric_u2_of_corresponding_continuous(
n2, a1, a2, d, c1, c2
):
return (a1 - a2) * d**2 * (n2 + c2) * (n2 + c2 + 1 / 2) - a1
[docs]def discrete_confocal_conics_ic(n1, n2, a1, a2, k, h, c1, c2):
"""discrete confocal conics -- parametrization with elliptic coordinates
leading to IC-nets"""
def A1(a1, a2, k, h):
return sqrt(a1 - a2) * dn(h / 2, k) / (k * cn(h / 2, k))
def f1(u1, u2, a1, a2, k, h, c1, c2):
return (
A1(a1, a2, k, h)
* (sn(h * n1 + c1, k) * dn(h * n2 + c2, k))
/ cn(h * n2 + c2, k)
)
def A2(a1, a2, k, h):
return sqrt(a1 - a2) * sqrt(1 - k**2) / (k * cn(h / 2, k))
def f2(u1, u2, a1, a2, k, h, c1, c2):
return A2(a1, a2, k, h) * cn(h * n1 + c1, k) / cn(h * n2 + c2, k)
return np.array(
[f1(n1, n2, a1, a2, k, h, c1, c2), f2(n1, n2, a1, a2, k, h, c1, c2)]
)
[docs]def discrete_confocal_conics_ic_u1_of_corresponding_continuous(
n1, a1, a2, k, h, c1, c2
):
return (a1 - a2) * (dn(h / 2, k) / cn(h / 2, k)) * sn(h * n1 + c1, k) * sn(
h * (n1 + 0.5) + c1, k
) - a1
[docs]def discrete_confocal_conics_ic_u2_of_corresponding_continuous(
n2, a1, a2, k, h, c1, c2
):
return (a1 - a2) * (dn(h / 2, k) / (k**2 * cn(h / 2, k))) * (
dn(h * n2 + c2, k) * dn(h * (n2 + 0.5) + c2, k)
) / (cn(h * n2 + c2, k) * cn(h * (n2 + 0.5) + c2, k)) - a1
[docs]def discrete_confocal_conics_ic_hyperbolic(n1, n2, a1, a2, k, h, c1, c2):
"""discrete confocal conics -- parametrization with elliptic coordinates
leading to IC-nets"""
def A1(a1, a2, k, h):
return k * sqrt(a1 - a2) * cn(h / 2, k) / dn(h / 2, k)
def f1(u1, u2, a1, a2, k, h, c1, c2):
return A1(a1, a2, k, h) * sn(h * n1 + c1, k) / sn(h * n2 + c2, k)
def A2(a1, a2, k, h):
return sqrt(a1 - a2) / dn(h / 2, k)
def f2(u1, u2, a1, a2, k, h, c1, c2):
return (
A2(a1, a2, k, h)
* dn(h * n1 + c1, k)
* cn(h * n2 + c2, k)
/ sn(h * n2 + c2, k)
)
return np.array(
[f1(n1, n2, a1, a2, k, h, c1, c2), f2(n1, n2, a1, a2, k, h, c1, c2)]
)
[docs]def discrete_confocal_conics_hyperbolic_pencil(n1, n2, d, c1, c2):
"""discrete confocal conics -- parametrization st diagonals are vertical
lines and dual wrt circles"""
def f1(n1, n2, d, c1, c2):
return exp(d * (n1 + n2 + c1 + c2))
def g1(n1, d, c1):
q = exp(-2 * d)
return sqrt(1 - q) * qgamma(-n1 + 0.5 - c1, q) / qgamma(-n1 - c1, q)
def g2(n2, d, c2):
q = exp(2 * d)
return sqrt(q - 1) * qgamma(n2 + 0.5 + c2, q) / qgamma(n2 + c2, q)
def f2(n1, n2, d, c1, c2):
return g1(n1, d, c1) * g2(n2, d, c2)
return np.array([f1(n1, n2, d, c1, c2), f2(n1, n2, d, c1, c2)])
[docs]def discrete_confocal_conics_hyperbolic_pencil_u1_of_corresponding_continuous(
n1, d, c1
):
return exp(d * (n1 + c1 + 0.25)) * exp(d * (n1 + 0.5 + c1 + 0.25)) - 2
[docs]def discrete_confocal_conics_hyperbolic_pencil_u2_of_corresponding_continuous(
n2, d, c2
):
return exp(d * (n2 + c2 - 0.25)) * exp(d * (n2 + 0.5 + c2 - 0.25)) - 2
[docs]def sampled_confocal_conic(
u, a1, a2, samples=10000, end=100, quadrants=None, prefer_ellipse=True
):
"""Coordinates of a sampled version of the conic
x^2 / (a1 + u) + y^2 / (a2 + u) = 1.
Returns two arrays. One for x coordinates and one for y coordinates.
"""
if u < -a1: # imaginary
return [[], []]
elif u < -a2 or (not prefer_ellipse and u == -a2): # hyperbola
if not quadrants:
u2_range = np.linspace(-a2 - end, -a2 + end, samples)
coords = [confocal_conics_sqrt(u, u2, a1, a2) for u2 in u2_range]
coords += [[np.nan, np.nan]]
coords += [confocal_conics_sqrt(-u - 2 * a1, u2, a1, a2) for u2 in u2_range]
else:
coords = []
for i in quadrants:
if i in {1, 2}:
u2_range = np.linspace(-a2, -a2 + end, samples)
elif i in {3, 4}:
u2_range = np.linspace(-a2 - end, -a2, samples)
if i in {1, 4}:
coords += [confocal_conics_sqrt(u, u2, a1, a2) for u2 in u2_range]
elif i in {2, 3}:
coords += [
confocal_conics_sqrt(-u - 2 * a1, u2, a1, a2) for u2 in u2_range
]
coords += [[np.nan, np.nan]]
return np.transpose(coords)
else: # ellipse
if not quadrants:
u1_range = np.linspace(-a1, -a1 + 4 * (a1 - a2), samples)
coords = [confocal_conics_sqrt(u1, u, a1, a2) for u1 in u1_range]
else:
coords = []
for i in quadrants:
if i == 4:
i = 2
elif i == 2:
i = 4
u1_range = np.linspace(
-a1 + (i - 1) * (a1 - a2), -a1 + i * (a1 - a2), samples
)
coords += [confocal_conics_sqrt(u1, u, a1, a2) for u1 in u1_range]
coords += [[np.nan, np.nan]]
return np.transpose(coords)