import numpy as np
from scipy.special import gamma
from ddg.datastructures.nets.domain import SmoothDomain, SmoothInterval
from ddg.datastructures.nets.net import NetCollection, SmoothCurve, SmoothNet
from ddg.geometry.quadrics import Pencil, Quadric
from ddg.math.functions import K, cd, cn, cs, dc, dn, ds, nc, nd, ns, sc, sd, sn
[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):
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)