# 思考：两个椭圆片能粘合成一个立体吗？

$$\frac{x^2}{a^2}+\frac{y^2}{b^2}=1$$

$$s=\int_0^y \sqrt{1+\left(\frac{dx}{dy}\right)^2}dy=\int_0^y \sqrt{\frac{a^2 y^2}{b^4-b^2 y^2}+1}dy\label{eq:bs}$$

$$y=h(s)$$

$$z=h(l-s)$$

$$\left(\frac{dx(s)}{ds}\right)^2+\left(\frac{dh(s)}{ds}\right)^2+\left(\frac{dh(l-s)}{ds}\right)^2=1$$

$$x(s) = \int_0^s \sqrt{1-\left(\frac{dh(s)}{ds}\right)^2-\left(\frac{dh(l-s)}{ds}\right)^2}ds\label{eq:f}$$

$$1-\left(\frac{dh(s)}{ds}\right)^2-\left(\frac{dh(l-s)}{ds}\right)^2 \geq 0\label{eq:c}$$

#! -*- coding: utf-8 -*-

import numpy as np

a = 2
b = 1
h = 1e-6 # 步长

def int(y, x):
"""自定义数值积分函数
"""
assert len(y) == len(x)
dx = np.diff(x)
_x = (x[:-1] + x[1:]) / 2
_y = (y[:-1] + y[1:]) / 2
s = np.cumsum(_y * dx)
s = np.interp(x, _x, s, left=0)
return s, x

def D(y, x):
"""自定义数值微分函数
"""
assert len(y) == len(x)
dy = np.diff(y) + 1e-10
dx = np.diff(x) + 1e-10
dy_dx = dy / dx
_x = (x[:-1] + x[1:]) / 2
dy_dx = np.interp(x, _x, dy_dx)
return dy_dx, x

def s2y():
"""弧长到椭圆y轴的映射（即h(s)。）
"""
y = np.arange(0, b, h)
ds = np.sqrt(1 + a**2 * y**2 / b**2 / (b**2 - y**2))
return int(ds, y)[::-1]

y, s = s2y()
l = s[-1]
z = np.interp(l - s, s, y)
dx2 = 1 - D(y, s)[0]**2 - D(z, s)[0]**2
assert (dx2 > - 10 * h).all() # 判断是否能粘合，即在10倍步长的误差内是否恒大于等于0
dx = np.sqrt(np.clip(dx2, 0, 1))
x = int(dx, s)[0]

# 采样1/10的点去绘图
x = x[::10]
y = y[::10]
z = z[::10]

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

plt.clf()
fig = plt.figure(figsize=(8,8))
ax = fig.gca(projection='3d')
ax.plot(x, y, z, color='blue')
ax.plot(x, - y, z, color='blue')
ax.plot(x, y, - z, color='green')
ax.plot(x, - y, - z, color='green')
ax.plot(x, np.zeros_like(x), z, color='red')
ax.plot(x, np.zeros_like(x), - z, color='red')

for i in range(0, len(x), 2000):
_ = ax.plot([x[i], x[i]], [y[i], - y[i]], [z[i], z[i]], color='blue')
_ = ax.plot([x[i], x[i]], [y[i], - y[i]], [- z[i], - z[i]], color='green')

for i in [-1, -1000]:
_ = ax.plot([x[i], x[i]], [y[i], - y[i]], [z[i], z[i]], color='blue')
if i != -1:
_ = ax.plot([x[i], x[i]], [y[i], - y[i]], [- z[i], - z[i]], color='green')

ax.set_xlim(0, a)
ax.set_ylim(-b * 1.2, b * 1.2)
ax.set_zlim(-b * 1.2, b * 1.2)

plt.show()