import numpy as np
from scipy.special import gamma
from ddg.math.functions import K, cd, cn, cs, dc, dn, ds, nc, nd, ns, sc, sd, sn
from ddg.nets._domain import SmoothDomain, SmoothInterval
from ddg.nets._net import NetCollection, SmoothCurve, SmoothNet
[docs]def discrete_confocal_coordinates_gamma(a=None, b=None, c=0):
if a is None:
if b is not None:
a = 2 * b
else:
a = 8
if b is None:
b = int(a / 2)
if c is None:
c = 0
intervals = [[-a, -b], [-b, -c], [-c, np.inf]]
def fct(n1, n2, n3, a=a, b=b, c=c):
A = 1 / (np.sqrt(a - b - 0.5) * np.sqrt(a - c - 1))
B = 1 / (np.sqrt(a - b - 0.5) * np.sqrt(b - c - 0.5))
C = 1 / (np.sqrt(a - c - 1) * np.sqrt(b - c - 0.5))
def f1(n1, n2, n3, a=a):
return (
A
* (gamma(n1 + a + 0.5) / gamma(n1 + a))
* (gamma(n2 + a) / gamma(n2 + a - 0.5))
* (gamma(n3 + a - 0.5) / gamma(n3 + a - 1.0))
)
def f2(n1, n2, n3, b=b):
return (
B
* (gamma(-n1 - b + 0.5) / gamma(-n1 - b))
* (gamma(n2 + b + 0.5) / gamma(n2 + b))
* (gamma(n3 + b) / gamma(n3 + b - 0.5))
)
def f3(n1, n2, n3, c=c):
return (
C
* (gamma(-n1 - c) / gamma(-n1 - c - 0.5))
* (gamma(-n2 - c + 0.5) / gamma(-n2 - c))
* (gamma(n3 + c + 0.5) / gamma(n3 + c))
)
return np.array((f1(n1, n2, n3), f2(n1, n2, n3), f3(n1, n2, n3)))
domain = SmoothDomain(intervals)
return SmoothNet(fct, domain)
[docs]def confocal_coordinates_elliptic(a=3, b=1, c=0):
k1 = np.sqrt((a - b) / (a - c))
k2 = np.sqrt(1 - k1**2)
k3 = k1
intervals = [[0.0, 4 * K(k1), True], [0.0, 4 * K(k2), True], [0.0, 2 * K(k3)]]
def fct(s1, s2, s3):
A = np.sqrt(a - c)
k1 = np.sqrt((a - b) / (a - c))
k2 = np.sqrt(1 - k1**2)
k3 = k1
return np.array(
(
A * sn(s1, k1) * dn(s2, k2) * ns(s3, k3),
A * cn(s1, k1) * cn(s2, k2) * ds(s3, k3),
A * dn(s1, k1) * sn(s2, k2) * cs(s3, k3),
)
)
domain = SmoothDomain(intervals)
return SmoothNet(fct, domain)
[docs]def confocal_coordinates_circles_normalized(a=None, b=None, c=None):
def sqrt(x):
if x < 0.0:
return 0.0
else:
return np.sqrt(x)
if a is None:
if b is not None:
a = 2 * b
else:
a = 8
if b is None:
b = int(a / 2)
if c is None:
c = 0
s0 = np.arccos(sqrt((a - b) / (a - c)))
intervals = [[s0, np.pi - s0], [-s0, s0], [0, np.pi / 2]]
def fct(s1, s2, s3, a=a, b=b, c=c):
A = (a - c) / (sqrt((a - b) * (b - c)))
B = 1 / (sqrt((a - b) * (b - c)))
C = (a - c) / (sqrt((a - b) * (b - c)))
def f1(s1, s2, s3):
return A * np.cos(s1) * np.cos(s2) * np.cos(s3)
def f2(s1, s2, s3, a=a, b=b, c=c):
return (
B
* sqrt(a - b - (a - c) * np.cos(s1) ** 2)
* sqrt((a - c) * np.cos(s2) ** 2 - a + b)
)
def f3(s1, s2, s3):
return C * np.sin(s1) * np.sin(s2) * np.sin(s3)
return np.array((f1(s1, s2, s3), f2(s1, s2, s3), f3(s1, s2, s3)))
domain = SmoothDomain(intervals)
return SmoothNet(fct, domain)
[docs]def minkowski_confocal_coordinates(k=None):
if k is None:
k = 0.8
KK = K(k)
intervals = [[0, 4 * KK, True], [-KK, KK], [0, 2 * KK]]
def fct(s1, s2, s3, k=k):
return np.array(
(
sn(s1, k) * dc(s2, k) * ns(s3, k),
cn(s1, k) * nc(s2, k) * ds(s3, k),
dn(s1, k) * sc(s2, k) * cs(s3, k),
)
)
domain = SmoothDomain(intervals)
net = SmoothNet(fct, domain)
return net
[docs]def minkowski_confocal_coordinates2(k=None):
if k is None:
k = 0.8
KK = K(k)
intervals = [[0, 4 * KK, True], [0, 4 * KK, True], [0, 2 * KK]]
def fct(s1, s2, s3, k=k):
return np.array(
(
sn(s1, k) * cd(s2, k) * ns(s3, k),
dn(s1, k) * nd(s2, k) * cs(s3, k) / k,
cn(s1, k) * sd(s2, k) * ds(s3, k),
)
)
domain = SmoothDomain(intervals)
return SmoothNet(fct, domain)
[docs]def minkowski_confocal_quadrics(a=1.0, b=0.75, c=0.0):
from ddg.geometry import Pencil, Quadric
Q1 = Quadric(np.diag([a, b, -c, -1.0]))
Q2 = Quadric(np.diag([1.0, 1.0, -1.0, 0.0]))
return Pencil(Q1, Q2)
[docs]def octahedral_grid_planes(k=None):
if k is None:
k = 0.5
KK = K(k)
interval = [0, 4 * KK, True]
def fct1(x):
return np.array((cn(x, k), -sn(x, k), 1, dn(x, k)), dtype="float")
def fct2(x):
return np.array((cn(x, k), -sn(x, k), 1, -dn(x, k)), dtype="float")
def fct3(x):
return np.array((-cn(x, k), sn(x, k), 1, -dn(x, k)), dtype="float")
def fct4(x):
return np.array((-cn(x, k), sn(x, k), 1, dn(x, k)), dtype="float")
domain = SmoothInterval(interval)
planes = [
SmoothCurve(fct1, domain),
SmoothCurve(fct2, domain),
SmoothCurve(fct3, domain),
SmoothCurve(fct4, domain),
]
return NetCollection(nets=planes)
[docs]def octahedral_grid_planes2(k=None):
if k is None:
k = 0.5
KK = K(k)
interval = [0, 4 * KK, True]
def fct1(x):
return np.array((dn(x, k), -k * sn(x, k), 1, cn(x, k)), dtype="float")
def fct2(x):
return np.array((dn(x, k), -k * sn(x, k), 1, -cn(x, k)), dtype="float")
def fct3(x):
return np.array((-dn(x, k), k * sn(x, k), 1, -cn(x, k)), dtype="float")
def fct4(x):
return np.array((-dn(x, k), k * sn(x, k), 1, cn(x, k)), dtype="float")
domain = SmoothInterval(interval)
planes = [
SmoothCurve(fct1, domain),
SmoothCurve(fct2, domain),
SmoothCurve(fct3, domain),
SmoothCurve(fct4, domain),
]
return NetCollection(nets=planes)