import numpy as np
[docs]def delta(N1, N2):
return np.pi / (2 * N1 + 2 * N2 - 1)
[docs]def q(N1, N2):
return np.cos((4 * N2 - 1) * np.pi / (4 * N1 + 4 * N2 - 2)) / np.cos(
np.pi / (4 * N1 + 4 * N2 - 2)
)
[docs]def abc(N1, N2, A):
c = 0.0
a = 2 * A / (q(N1, N2) + 1)
b = (1 - q(N1, N2)) * a / 2
return a, b, c
[docs]def g1_function(a, b, c, delta, g1_0):
def f1(n1):
return np.sqrt((a - c) / np.cos(delta / 2)) * np.sin(delta * n1)
def g1(n1):
if n1 == 0:
return g1_0
elif n1 > 0:
return (a - b - f1(n1 - 0.5) * f1(n1)) / g1(n1 - 0.5)
elif n1 < 0:
return (a - b - f1(n1) * f1(n1 + 0.5)) / g1(n1 + 0.5)
return g1
[docs]def g2_function(a, b, c, delta, g2_0):
def f2(n2):
return np.sqrt((a - c) / np.cos(delta / 2)) * np.cos(delta * n2)
def g2(n2):
if n2 == 0:
return g2_0
elif n2 > 0:
return (f2(n2 - 0.5) * f2(n2) - a + b) / g2(n2 - 0.5)
elif n2 < 0:
return (f2(n2) * f2(n2 + 0.5) - a + b) / g2(n2 + 0.5)
return g2
[docs]def ellipsoid(n1, n2, s3, N1, N2, A):
d = delta(N1, N2)
eps = 1 / np.cos(d / 2)
a, b, c = abc(N1, N2, A)
sqrt = np.sqrt((a - b) * (b - c))
g1_0 = np.sqrt(a - b)
g2_0 = np.sqrt(b - c)
if n1 <= N1:
g1 = g1_function(a, b, c, d, g1_0)
else:
g1 = g1_function(a, b, c, d, -g1_0)
n1 = 2 * N1 - n1
g2 = g2_function(a, b, c, d, g2_0)
x = eps * np.sin(d * n1) * np.cos(d * n2) * np.cos(s3) * (a - c) / sqrt
y = g1(n1) * g2(n2) / sqrt
z = eps * np.cos(d * n1) * np.sin(d * n2) * np.sin(s3) * (a - c) / sqrt
return np.array([x, y, z])