"""Discrete parametrized surfaces and coordinate systems."""
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
[docs]def discrete_confocal_conics_gamma(n1, n2, b1, b2):
r"""2-dimensional discrete confocal coordinate system defined on `1/2 Z^2`.
The parametrization is given in terms of gamma functions ::
x(n1, n2) = x_sgn * D * (n1 + b1)_1/2 * (n2 + b1 - 0.5)_1/2
y(n1, n2) = y_sgn * D * (-n1 - b2)_1/2 * (n2 + b2)_1/2
where ::
D = 1 / sqrt(b1 - b2 - 0.5)
and (.)_1/2 denotes the Pochhammer symbol ::
(n)_1/2 = gamma(n + 1/2) / gamma(n).
Incident points from dual sublattices, e.g. `Z^2` and `(Z+1/2)^2`
are related by polarity ::
x * x* / (a1 + u) + y * y* / (a2 + u) = 1
where ::
a1 = b1 + 0.5
a2 = b2 + 1.
(d1, d2) in {-1,1} x {-1, 1}
x = x(n1, n2)
x* = x(n1 + 0.5*d1, n2 + 0.5*d2)
y = y(n1, n2)
y* = y(n1 + 0.5*d1, n2 + 0.5*d2)
u = u1(n1 + 0.25*d1) or u = u2(n2 + 0.25*d2)
The relation to the corresponding confocal conics that constitute the polarity is
given by ::
u1(n1 + 0.25) = n1 + b1 - a1 = n1 - 0.5
u2(n2 + 0.25) = n2 + b2 - a2 = n2 - 1.0
as computed in `discrete_confocal_conics_gamma_u1`
and `discrete_confocal_conics_gamma_u2`.
With `x_sgn = y_sgn = 1` this defines a parametrization of the first quadrant.
Here the parametrization is extended to the entire plane by reflection along the
coordinate axes. This is done by computing `x_sgn` and `y_sgn` accordingly.
Parameters
----------
n1 : float in [-b1, -b2]
First parameter (should be integer and half-integer values).
n2 : float in [-b2, inf)
Second parameter (should be integer and half-integer values).
b1 : float with b1 > b2
Determines the first semi axis of the conics.
b2 : float with b1 > b2
Determines the second semi axis of the conics.
Returns
-------
np.ndarray of shape (2,)
Point in R^2.
"""
x1_sign = 1
x2_sign = 1
n1 = ((n1 + b1) % (4 * (b1 - b2))) - b1
if -b2 <= n1 < b1 - 2 * b2:
n1 = -n1 - 2 * b2
x2_sign *= -1
elif b1 - 2 * b2 <= n1 < 2 * b1 - 3 * b2:
n1 = n1 - 2 * (b1 - b2)
x1_sign *= -1
x2_sign *= -1
elif 2 * b1 - 3 * b2 <= n1 < 3 * b1 - 4 * b2:
n1 = -n1 + 2 * b1 - 4 * b2
x1_sign *= -1
if n2 < -b2:
n2 = -n2 - 2 * b2
x2_sign *= -1
def D(b1, b2):
return 1 / sqrt(b1 - b2 - 0.5)
def f1(n1, n2, b1, b2):
return (
x1_sign
* D(b1, b2)
* (gamma(n1 + b1 + 0.5) / gamma(n1 + b1))
* (gamma(n2 + b1) / gamma(n2 + b1 - 0.5))
)
def f2(n1, n2, b1, b2):
return (
x2_sign
* D(b1, b2)
* (gamma(-n1 - b2 + 0.5) / gamma(-n1 - b2))
* (gamma(n2 + b2 + 0.5) / gamma(n2 + b2))
)
return np.array([f1(n1, n2, b1, b2), f2(n1, n2, b1, b2)])
[docs]def discrete_confocal_conics_gamma_u1(n1):
r"""Auxiliary function for discrete confocal coordinates
`discrete_confocal_conics_gamma` to compute the cooresponding continuous
confocal conics ::
x**2 / (a1 + u) + y**2 / (a2 + u) = 1
that establish the polarity relation.
It computes the `u = u1` value for the conic adjacent to points
with parameter `n1`: ::
u1(n1 + 0.25) = n1 - 0.5
Note that `0.25` is implicitely added to the input parameter.
Parameters
----------
n1 : float
Should be integer or half-integer value.
Returns
-------
float
The corresponding u1 value.
See Also
--------
discrete_confocal_conics_gamma
"""
return n1 - 0.5
[docs]def discrete_confocal_conics_gamma_u2(n2):
r"""Auxiliary function for discrete confocal coordinates
`discrete_confocal_conics_gamma` to compute the cooresponding continuous
confocal conics ::
x**2 / (a1 + u) + y**2 / (a2 + u) = 1
that establish the polarity relation.
It computes the `u = u2` value for the conic adjacent to points
with parameter `n2`: ::
u2(n2 + 0.25) = n2 - 1.0
Note that `0.25` is implicitely added to the input parameter.
Parameters
----------
n2 : float
Should be integer or half-integer value.
Returns
-------
float
The corresponding u2 value.
See Also
--------
discrete_confocal_conics_gamma
"""
return n2 - 1.0
[docs]def discrete_confocal_conics_trigonometric(n1, n2, a1, a2, h1, h2, c1=0.0, c2=0.0):
"""2-dimensional discrete confocal coordinate system defined on `1/2 Z^2`.
The parametrization is given in terms of trigonometric functions ::
x(n1, n2) = A * cos(h1 * n1 + c1) * cosh(h2 * n2 + c2)
y(n1, n2) = A * sin(h1 * n1 + c1) * sinh(h2 * n2 + c2)
where ::
A = sqrt((a1 - a2) / (cos(h1 / 2) * cosh(h2 / 2))).
Incident points from dual sublattices are related by polarity.
The parameters of the corresponding confocal conics that constitute the polarity are
computed in `discrete_confocal_conics_trigonometric_u1`
and `discrete_confocal_conics_trigonometric_u2`.
The parametrization becomes periodic in `n1` for `h1 = 2*pi / m`
with some integer `m`.
Parameters
----------
n1 : float
First parameter (should be integer and half-integer values).
n2 : float
Second parameter (should be integer and half-integer values).
a1 : float with a1 > a2
Determines the first semi axis of the conics.
a2 : float with a1 > a2
Determines the second semi axis of the conics.
h1 : float
Frequency in first parameter direction.
h2 : float
Frequency in second parameter direction.
c1 : float
Shift in first parameter direction.
c2 : float
Shift in second parameter direction.
Returns
-------
np.ndarray of shape (2,)
Point in R^2.
"""
def A(a1, a2, h1, h2):
return sqrt((a1 - a2) / (cos(h1 / 2) * cosh(h2 / 2)))
def x1(n1, n2, a1, a2, h1, h2, c1, c2):
return A(a1, a2, h1, h2) * cos(h1 * n1 + c1) * cosh(h2 * n2 + c2)
def x2(n1, n2, a1, a2, h1, h2, c1, c2):
return A(a1, a2, h1, h2) * sin(h1 * n1 + c1) * sinh(h2 * n2 + c2)
return np.array(
[x1(n1, n2, a1, a2, h1, h2, c1, c2), x2(n1, n2, a1, a2, h1, h2, c1, c2)]
)
[docs]def discrete_confocal_conics_trigonometric_u1(n1, a1, a2, h1, c1=0.0):
"""Auxiliary function for discrete confocal coordinates
`discrete_confocal_conics_trigonometric` to compute the cooresponding continuous
confocal conics that establish the polarity relation.
It computes the `u = u1(n1 + 0.25)` value for the conic adjacent to points
with parameter `n1`.
Parameters
----------
n1 : float
Should be integer or half-integer value.
a1 : float with a1 > a2
Determines the first semi axis of the conics.
a2 : float with a1 > a2
Determines the second semi axis of the conics.
h1 : float
Frequency in first parameter direction.
c1 : float
Shift in first parameter direction.
Returns
-------
float
The corresponding u1 value.
See Also
--------
discrete_confocal_conics_trigonomitric
"""
return (
-a1 * sin(h1 * n1 + c1) ** 2
- a2 * cos(h1 * n1 + c1) ** 2
- (a1 - a2) * sin(h1 * n1 + c1) * cos(h1 * n1 + c1) * tan(h1 / 2)
)
[docs]def discrete_confocal_conics_trigonometric_u2(n2, a1, a2, h2, c2=0.0):
"""Auxiliary function for discrete confocal coordinates
`discrete_confocal_conics_trigonometric` to compute the cooresponding continuous
confocal conics that establish the polarity relation.
It computes the `u = u2(n2 + 0.25)` value for the conic adjacent to points
with parameter `n2`.
Parameters
----------
n2 : float
Should be integer or half-integer value.
a1 : float with a1 > a2
Determines the first semi axis of the conics.
a2 : float with a1 > a2
Determines the second semi axis of the conics.
h2 : float
Frequency in second parameter direction.
c2 : float
Shift in second parameter direction.
Returns
-------
float
The corresponding u2 value.
See Also
--------
discrete_confocal_conics_trigonomitric
"""
return (
a1 * sinh(h2 * n2 + c2) ** 2
- a2 * cosh(h2 * n2 + c2) ** 2
+ (a1 - a2) * sinh(h2 * n2 + c2) * cosh(h2 * n2 + c2) * tanh(h2 / 2)
)
[docs]def discrete_confocal_conics_concentric(n1, n2, a1, a2, d, c1=0.0, c2=0.0):
"""2-dimensional discrete confocal coordinate system defined on `1/2 Z^2`.
The parametrization is given by ::
x(n1, n2) = sqrt(a1 - a2) * d**2 * (n1 + c1) * (n2 + c2)
y(n1, n2) = sqrt(a1 - a2) * (n1 + p(c1))_1/2 * (-(n1 + m(c1)) + 0.5)_1/2
* (n2 + p(c2))_1/2 * (-(n2 + m(c2)))_1/2
where (.)_1/2 denotes the Pochhammer symbol ::
(n)_1/2 = gamma(n + 1/2) / gamma(n),
and ::
p(c) = c + 0.25 + sqrt(16 + d**2) / (4*d),
m(c) = c + 0.25 - sqrt(16 + d**2) / (4*d).
Incident points from dual sublattices are related by polarity.
The parameters of the corresponding confocal conics that constitute the polarity are
computed in `discrete_confocal_conics_concentric_u1`
and `discrete_confocal_conics_concentric_u2`.
With c1 = c2 = 0 and ::
d = 4 / sqrt( (2*l + 1)**2 - 1 ) for some positive integer l
one can achieve boundary conditions for `n1` in `[-(l+1)/2, (l+1)/2]`
and `n2` in `[l/2, inf]`.
Parameters
----------
n1 : float
First parameter (should be integer and half-integer values).
n2 : float
Second parameter (should be integer and half-integer values).
a1 : float with a1 > a2
Determines the first semi axis of the conics.
a2 : float with a1 > a2
Determines the second semi axis of the conics.
d : float
Scaling factor.
c1 : float
Shift in first parameter direction.
c2 : float
Shift in second parameter direction.
Returns
-------
np.ndarray of shape (2,)
Point in R^2.
Notes
-----
Not tested yet.
"""
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 (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 (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) * d**2 * g1(n1, d, c1) * g2(n2, d, c2)
return np.array([x1(n1, n2, a1, a2, d, c1, c2), x2(n1, n2, a1, a2, d, c1, c2)])
[docs]def discrete_confocal_conics_concentric_u1(n1, a1, a2, d, c1=0):
"""Auxiliary function for discrete confocal coordinates
`discrete_confocal_conics_concentric` to compute the cooresponding continuous
confocal conics that establish the polarity relation.
Parameters
----------
n1 : float
a1 : float
a2 : float
d : float
c1 : float
Returns
-------
float
The corresponding u1 value.
See Also
--------
discrete_confocal_conics_concentric
"""
return (a1 - a2) * d**2 * (n1 + c1) * (n1 + c1 + 1 / 2) - a1
[docs]def discrete_confocal_conics_concentric_u2(n2, a1, a2, d, c2=0):
"""Auxiliary function for discrete confocal coordinates
`discrete_confocal_conics_concentric` to compute the cooresponding continuous
confocal conics that establish the polarity relation.
Parameters
----------
n1 : float
a1 : float
a2 : float
d : float
c2 : float
Returns
-------
float
The corresponding u2 value.
See Also
--------
discrete_confocal_conics_concentric
"""
return (a1 - a2) * d**2 * (n2 + c2) * (n2 + c2 + 1 / 2) - a1
[docs]def discrete_confocal_conics_ic_ellipse(n1, n2, a1, a2, k, h, c1=0.0, c2=0.0):
"""2-dimensional discrete confocal coordinate system defined on `1/2 Z^2`.
The parametrization is given by ::
x(n1, n2) = A1 * (sn(h * n1 + c1, k) * dn(h * n2 + c2, k)) / cn(h * n2 + c2, k)
y(n1, n2) = A2 * cn(h * n1 + c1, k) / cn(h * n2 + c2, k)
where ::
A1 = sqrt(a1 - a2) * dn(h / 2, k) / (k * cn(h / 2, k))
A2 = sqrt(a1 - a2) * sqrt(1 - k**2) / (k * cn(h / 2, k))
Incident points from dual sublattices are related by polarity.
The parameters of the corresponding confocal conics that constitute the polarity are
computed in `discrete_confocal_conics_ic_ellipse_u1`
and `discrete_confocal_conics_ic_ellipse_u2`.
The constants `k` and `h` should satisfy `0 < k**2 < 1` and `cn(h/2,k) > 0`.
Symmetry is achieved by `c1 = c1 = 0` and periodicity by ::
h = K(k) / m for some positive integer or half-integer m.
Then the parameters may be restricted to ::
n1 in [-2*m, 2*m], n2 in [0, m - 0.5].
Parameters
----------
n1 : float
First parameter (should be integer and half-integer values).
n2 : float
Second parameter (should be integer and half-integer values).
a1 : float with a1 > a2
Determines the first semi axis of the conics.
a2 : float with a1 > a2
Determines the second semi axis of the conics.
k : float
Modulus of the elliptic functions (0 < k**2 < 1).
h : float
Scaling / frequency factor
c1 : float
Shift in first parameter direction.
c2 : float
Shift in second parameter direction.
Returns
-------
np.ndarray of shape (2,)
Point in R^2.
Notes
-----
Not tested yet.
"""
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_ellipse_u1(n1, a1, a2, k, h, c1):
"""Auxiliary function for discrete confocal coordinates
`discrete_confocal_conics_ic_ellipse` to compute the cooresponding continuous
confocal conics that establish the polarity relation.
Parameters
----------
n1 : float
a1 : float
a2 : float
k : float
h : float
c1 : float
Returns
-------
float
The corresponding u1 value.
See Also
--------
discrete_confocal_conics_ic_ellipse
"""
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_ellipse_u2(n2, a1, a2, k, h, c2):
"""Auxiliary function for discrete confocal coordinates
`discrete_confocal_conics_ic_ellipse` to compute the cooresponding continuous
confocal conics that establish the polarity relation.
Parameters
----------
n1 : float
a1 : float
a2 : float
k : float
h : float
c2 : float
Returns
-------
float
The corresponding u2 value.
See Also
--------
discrete_confocal_conics_ic_ellipse
"""
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_hyperbola(n1, n2, a1, a2, k, h, c1, c2):
"""2-dimensional discrete confocal coordinate system defined on `1/2 Z^2`.
The parametrization is given by ::
x(n1, n2) = A1 * sn(h * n1 + c1, k) / sn(h * n2 + c2, k)
y(n1, n2) = A2 * dn(h * n1 + c1, k) * cn(h * n2 + c2, k) / sn(h * n2 + c2, k)
where ::
A1 = k * sqrt(a1 - a2) * cn(h / 2, k) / dn(h / 2, k)
A2 = sqrt(a1 - a2) / dn(h / 2, k)
Incident points from dual sublattices are related by polarity.
Parameters
----------
n1 : float
First parameter (should be integer and half-integer values).
n2 : float
Second parameter (should be integer and half-integer values).
a1 : float with a1 > a2
Determines the first semi axis of the conics.
a2 : float with a1 > a2
Determines the second semi axis of the conics.
k : float
Modulus of the elliptic functions (0 < k**2 < 1).
h : float
Scaling / frequency factor
c1 : float
Shift in first parameter direction.
c2 : float
Shift in second parameter direction.
Returns
-------
np.ndarray of shape (2,)
Point in R^2.
Notes
-----
Not tested yet.
"""
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=0.0, c2=0.0):
"""2-dimensional discrete confocal coordinate system defined on `1/2 Z^2`.
The parametrization is given by ::
x(n1, n2) = exp(d * (n1 + n2 + c1 + c2))
y(n1, n2) = sqrt(1 - 1/q) * (qgamma(-n1 + 0.5 - c1, 1/q) / qgamma(-n1 -c1, 1/q))
* sqrt(q - 1) * qgamma(n2 + 0.5 + c2, q) / qgamma(n2 + c2, q)
where ::
q = exp(2*d)
and `qgamma` denotes the q-gamma function, whiche is characerized by
qgamme(z+1, q) = qgamma(z, q) * (1 - q**z) / (1 - q)
Incident points from dual sublattices are related by polarity.
The parameters of the corresponding confocal conics that constitute the polarity are
computed in `discrete_confocal_conics_hyperbolic_pencil_u1`
and `discrete_confocal_conics_ic_hyperbolic_pencil_u2`.
Boundary conditions may be achieved by `c1 = c2 = 0`,
and then the parameters may be restricted to `n1 <= 0`, `n2 >=0`.
Parameters
----------
n1 : float
First parameter (should be integer and half-integer values).
n2 : float
Second parameter (should be integer and half-integer values).
d : float
Scaling / frequency factor
c1 : float
Shift in first parameter direction.
c2 : float
Shift in second parameter direction.
Returns
-------
np.ndarray of shape (2,)
Point in R^2.
Notes
-----
Not tested yet.
"""
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(n1, d, c1):
"""Auxiliary function for discrete confocal coordinates
`discrete_confocal_conics_hyperbolic_pencil` to compute the cooresponding continuous
confocal conics that establish the polarity relation.
Parameters
----------
n1 : float
d : float
c1 : float
Returns
-------
float
The corresponding u1 value.
See Also
--------
discrete_confocal_conics_ic_hyperbolic_pencil
"""
return exp(d * (n1 + c1 + 0.25)) * exp(d * (n1 + 0.5 + c1 + 0.25)) - 2
[docs]def discrete_confocal_conics_hyperbolic_pencil_u2(n2, d, c2):
"""Auxiliary function for discrete confocal coordinates
`discrete_confocal_conics_hyperbolic_pencil` to compute the cooresponding continuous
confocal conics that establish the polarity relation.
Parameters
----------
n2 : float
d : float
c2 : float
Returns
-------
float
The corresponding u2 value.
See Also
--------
discrete_confocal_conics_ic_hyperbolic_pencil
"""
return exp(d * (n2 + c2 - 0.25)) * exp(d * (n2 + 0.5 + c2 - 0.25)) - 2