Ответ на вопрос https://ru.stackoverflow.com/questions/1293488/

Пусть дана точка $P(x_p,y_p,z_p)$ и эллипсоид, у которого известны полуоси и центр, его уравнение: $x^2/a^2+y^2/a^2+z^2/b^2=1$. 

Точка $P$ лежит за пределами эллипсоида в некотором отдалении. Множество касательных к эллипсоиду, проведенных из $P$ образует коническую поверхность. Надо найти общий вид точек касания этого конуса и эллипсоида (понятно, что эти точки будут принадлежать эллипсу, то есть задача по сути сводится к отысканию уравнения эллипса касания конуса и эллипсоида). Можно ли здесь упростить задачу, перейдя к сферическим координатам? Если да, то как действовать дальше? Если нет, то вопрос аналогичный.


```python
import sympy as sp
import sympy.abc as abc
import sympy.vector as spv
```

$a$ и $b$ - параметры эллипсоида $\frac{x^2}{a^2} + \frac{y^2}{a^2} + \frac{z^2}{b^2}$


```python
a,b = abc.a, abc.b

```

Три системы координат:
- $A$ - исходная система координат, в которой задан эллипсоид
- $B$ - система координат, в которой эллипсоид превращается в единичную сферу
- $C$ - система координат, в которой эллипсоид превращен в сферу, а точка $P$ находится на оси $Z$

Матрицы:
- $A2B$ переводит координаты из $A$ в $B$, а матрица $B2A$ переводит обратно
- $B2C$ переводит координаты из $B$ в $C$, а матрица $C2B$ переводит обратно



```python
A2B = sp.Matrix([[1/a, 0,   0],
                 [0,   1/a, 0],
                 [0,   0,   1/b]])
```


```python
B2A = sp.Matrix([[a, 0, 0],
                 [0, a, 0],
                 [0, 0, b]])
```

Проверка того, что матрицы $A2B$ и $B2A$ обратны


```python
B2A*A2B
```




$\displaystyle \left[\begin{matrix}1 & 0 & 0\\0 & 1 & 0\\0 & 0 & 1\end{matrix}\right]$



В системе координат $B$ вектор $OP$ составляет с осью $x$ угол $\phi$, а с осью $\z$ угол $\theta$


Длина вектора $OP$ в этой системе координат равно $R = \sqrt{\frac{x^2}{a^2} + \frac{y^2}{a^2} + \frac{z^2}{b^2}}$, где $x,y,z$ - координаты точки $P$ в исходной системе координат.

Углы $\theta$ и $\phi$ в системе $B$ определяются так:
- $\cos\theta = \frac{z}{bR}$
- $\cos\phi = \frac{x}{aRsin\theta}$
- $\sin\phi = \frac{y}{aRsin\theta}$


```python
phi, theta = abc.phi, abc.theta
```

Поворот системы координат от $B$ к $C$, в которой точка $P$ окажется на оси $Z$, выполняется в два приёма:
- сначала на угол $\phi$ вокруг оси $Z$. В результате точка $P$ окажется в плоскости $ZOX$
- затем на угол $\theta$ вокруг оси $Y$. В результате точка $P$ окажется на оси $Z$ в новой системе координат.


```python
B2C = spv.BodyOrienter(phi, theta, 0, 'ZYZ').rotation_matrix()
B2C
```




$\displaystyle \left[\begin{matrix}\cos{\left(\phi \right)} \cos{\left(\theta \right)} & \sin{\left(\phi \right)} \cos{\left(\theta \right)} & - \sin{\left(\theta \right)}\\- \sin{\left(\phi \right)} & \cos{\left(\phi \right)} & 0\\\sin{\left(\theta \right)} \cos{\left(\phi \right)} & \sin{\left(\phi \right)} \sin{\left(\theta \right)} & \cos{\left(\theta \right)}\end{matrix}\right]$



Матрица обратного преобразования строится в обратном порядке: 
- сначала на угол $-\theta$ вокруг оси $Y$,
- затем на угол $-\phi$ вокруг оси $Z$. 


```python
C2B = spv.BodyOrienter(-theta, -phi, 0, 'YZX').rotation_matrix()
C2B
```




$\displaystyle \left[\begin{matrix}\cos{\left(\phi \right)} \cos{\left(\theta \right)} & - \sin{\left(\phi \right)} & \sin{\left(\theta \right)} \cos{\left(\phi \right)}\\\sin{\left(\phi \right)} \cos{\left(\theta \right)} & \cos{\left(\phi \right)} & \sin{\left(\phi \right)} \sin{\left(\theta \right)}\\- \sin{\left(\theta \right)} & 0 & \cos{\left(\theta \right)}\end{matrix}\right]$



Проверка того, что матрицы $B2C$ и $C2B$ обратны друг другу


```python
sp.trigsimp(C2B*B2C) # упрощение тригонометрических выражений
```




$\displaystyle \left[\begin{matrix}1 & 0 & 0\\0 & 1 & 0\\0 & 0 & 1\end{matrix}\right]$




```python
R = abc.R
```

Точка $P$ в системе координат $C$ находится на оси $Z$ на расстоянии $R$ от начала координат


```python
P_C = sp.Matrix([[0],[0],[R]])
```

Координаты точки $P$ в исходной системе координат


```python
P_A = B2A*C2B*P_C
```


```python
P_A
```




$\displaystyle \left[\begin{matrix}R a \sin{\left(\theta \right)} \cos{\left(\phi \right)}\\R a \sin{\left(\phi \right)} \sin{\left(\theta \right)}\\R b \cos{\left(\theta \right)}\end{matrix}\right]$



В системе координат $C$ эллипс представлен единичной сферой. Касательные из точки $P$ в этой системе координат образуют круг радиусом $r$ в плоскости, отстоящей на расстояние $h$ от начала координат.

Подробнее о касательной к кругу можно прочитать здесь: https://www.wikiwand.com/en/Tangent_lines_to_circles


```python
r = sp.sqrt(1 - 1/(R*R))
r
```




$\displaystyle \sqrt{1 - \frac{1}{R^{2}}}$




```python
h = sp.sqrt(1 - r*r)
```


```python
h
```




$\displaystyle \sqrt{\frac{1}{R^{2}}}$




```python

```

Построим уравнение касательного круга в параметрическом виде, как вращение точки вокруг оси $Z$ с круговой частотой $\omega$


```python
omega = abc.omega
t = abc.t
```


```python
base = sp.Matrix([[0],[0], [1/R]]) # центр круга
e1 = sp.Matrix([[r],[0],[0]]) # проекция точки на ось x
e2 = sp.Matrix([[0], [r],[0]]) # проекция точки на ось y
```

Радиус-вектор точек, в которых конус касается эллипсоида, в системе координат $C$


```python
tangent = base + e1*sp.cos(t) + e2*sp.sin(t)
tangent
```




$\displaystyle \left[\begin{matrix}\sqrt{1 - \frac{1}{R^{2}}} \cos{\left(t \right)}\\\sqrt{1 - \frac{1}{R^{2}}} \sin{\left(t \right)}\\\frac{1}{R}\end{matrix}\right]$



Проверка: касательная плоскость к единичной сфере в точке $(x_0, y_0, z_0)$ задаётся уравнением $2x_0(x-x_0) + 2y_0(y-y_0) + 2z_0(z-z_0) = 0$

Подставим в это уравнение в качестве $(x_0, y_0, z_0)$ координаты из $tangent$, а вместо $(x, y, z)$ координаты точки $P = (0, 0, R)$

Для справки уравнение касательной плоскости к поверхности, заданной уравнением $F(x,y,z)=0$ в точке $(x_0, y_0, z_0)$:

$F'_x(x_0)\cdot(x-x_0) + F'_y(y_0)\cdot(y-y_0) + F'_z(z_0)\cdot(z-z_0) = 0$


```python
_t = tangent
_u = _t[0]*(-_t[0]) + _t[1]*(-_t[1]) + _t[2]*(R-_t[2])
sp.trigsimp(_u)
```




$\displaystyle 0$



Радиус вектор точек , в которых конус касается эллипсоида, в системе координат $A$


```python
tangent_A = B2A*C2B*tangent
```


```python
tangent_A
```




$\displaystyle \left[\begin{matrix}- a \sqrt{1 - \frac{1}{R^{2}}} \sin{\left(\phi \right)} \sin{\left(t \right)} + a \sqrt{1 - \frac{1}{R^{2}}} \cos{\left(\phi \right)} \cos{\left(t \right)} \cos{\left(\theta \right)} + \frac{a \sin{\left(\theta \right)} \cos{\left(\phi \right)}}{R}\\a \sqrt{1 - \frac{1}{R^{2}}} \sin{\left(\phi \right)} \cos{\left(t \right)} \cos{\left(\theta \right)} + a \sqrt{1 - \frac{1}{R^{2}}} \sin{\left(t \right)} \cos{\left(\phi \right)} + \frac{a \sin{\left(\phi \right)} \sin{\left(\theta \right)}}{R}\\- b \sqrt{1 - \frac{1}{R^{2}}} \sin{\left(\theta \right)} \cos{\left(t \right)} + \frac{b \cos{\left(\theta \right)}}{R}\end{matrix}\right]$



Проверим, что найденные точки лежат на эллипсоиде


```python
_t = tangent_A
_u = _t[0]**2/a**2 + _t[1]**2/a**2 + _t[2]**2/b**2
sp.trigsimp(_u)
```




$\displaystyle 1$



Проверка 2: касательная плоскость к эллипсоиду в точке $(x_0, y_0, z_0)$ задаётся уравнением $2x_0/a^2(x-x_0) + 2y_0/a^2(y-y_0) + 2z_0/b^2(z-z_0) = 0$


```python
_t = tangent_A
_p = P_A
_u = _t[0]/a**2*(_p[0]-_t[0]) + _t[1]/a**2*(_p[1]-_t[1]) + _t[2]/b**2*(_p[2]-_t[2])
sp.trigsimp(_u)
```




$\displaystyle 0$



Подставим формулы для $R$ и углов $\theta$ и $\phi$ в выражение для точек касания конуса к эллипсу


```python
x,y,z = sp.symbols("x y z")
```


```python
t_a = tangent_A.subs({                 
                sp.cos(phi) : x/a/sp.sin(theta)/R,
                sp.sin(phi) : y/a/sp.sin(theta)/R
               })                                        \
    .subs({sp.sin(theta) : sp.sqrt(1-sp.cos(theta)**2)}) \
    .subs({sp.cos(theta) : z/b/R})                       \
    .subs({R:sp.sqrt((x/a)**2 + (y/a)**2 + (z/b)**2)})
t_a
```




$\displaystyle \left[\begin{matrix}\frac{x}{\frac{z^{2}}{b^{2}} + \frac{x^{2}}{a^{2}} + \frac{y^{2}}{a^{2}}} - \frac{y \sqrt{1 - \frac{1}{\frac{z^{2}}{b^{2}} + \frac{x^{2}}{a^{2}} + \frac{y^{2}}{a^{2}}}} \sin{\left(t \right)}}{\sqrt{1 - \frac{z^{2}}{b^{2} \left(\frac{z^{2}}{b^{2}} + \frac{x^{2}}{a^{2}} + \frac{y^{2}}{a^{2}}\right)}} \sqrt{\frac{z^{2}}{b^{2}} + \frac{x^{2}}{a^{2}} + \frac{y^{2}}{a^{2}}}} + \frac{x z \sqrt{1 - \frac{1}{\frac{z^{2}}{b^{2}} + \frac{x^{2}}{a^{2}} + \frac{y^{2}}{a^{2}}}} \cos{\left(t \right)}}{b \sqrt{1 - \frac{z^{2}}{b^{2} \left(\frac{z^{2}}{b^{2}} + \frac{x^{2}}{a^{2}} + \frac{y^{2}}{a^{2}}\right)}} \left(\frac{z^{2}}{b^{2}} + \frac{x^{2}}{a^{2}} + \frac{y^{2}}{a^{2}}\right)}\\\frac{x \sqrt{1 - \frac{1}{\frac{z^{2}}{b^{2}} + \frac{x^{2}}{a^{2}} + \frac{y^{2}}{a^{2}}}} \sin{\left(t \right)}}{\sqrt{1 - \frac{z^{2}}{b^{2} \left(\frac{z^{2}}{b^{2}} + \frac{x^{2}}{a^{2}} + \frac{y^{2}}{a^{2}}\right)}} \sqrt{\frac{z^{2}}{b^{2}} + \frac{x^{2}}{a^{2}} + \frac{y^{2}}{a^{2}}}} + \frac{y}{\frac{z^{2}}{b^{2}} + \frac{x^{2}}{a^{2}} + \frac{y^{2}}{a^{2}}} + \frac{y z \sqrt{1 - \frac{1}{\frac{z^{2}}{b^{2}} + \frac{x^{2}}{a^{2}} + \frac{y^{2}}{a^{2}}}} \cos{\left(t \right)}}{b \sqrt{1 - \frac{z^{2}}{b^{2} \left(\frac{z^{2}}{b^{2}} + \frac{x^{2}}{a^{2}} + \frac{y^{2}}{a^{2}}\right)}} \left(\frac{z^{2}}{b^{2}} + \frac{x^{2}}{a^{2}} + \frac{y^{2}}{a^{2}}\right)}\\- b \sqrt{1 - \frac{z^{2}}{b^{2} \left(\frac{z^{2}}{b^{2}} + \frac{x^{2}}{a^{2}} + \frac{y^{2}}{a^{2}}\right)}} \sqrt{1 - \frac{1}{\frac{z^{2}}{b^{2}} + \frac{x^{2}}{a^{2}} + \frac{y^{2}}{a^{2}}}} \cos{\left(t \right)} + \frac{z}{\frac{z^{2}}{b^{2}} + \frac{x^{2}}{a^{2}} + \frac{y^{2}}{a^{2}}}\end{matrix}\right]$




```python
center = B2A*C2B*base
v1 = B2A*C2B*e1
v2 = B2A*C2B*e2
```


```python
center
```




$\displaystyle \left[\begin{matrix}\frac{a \sin{\left(\theta \right)} \cos{\left(\phi \right)}}{R}\\\frac{a \sin{\left(\phi \right)} \sin{\left(\theta \right)}}{R}\\\frac{b \cos{\left(\theta \right)}}{R}\end{matrix}\right]$




```python
v1
```




$\displaystyle \left[\begin{matrix}a \sqrt{1 - \frac{1}{R^{2}}} \cos{\left(\phi \right)} \cos{\left(\theta \right)}\\a \sqrt{1 - \frac{1}{R^{2}}} \sin{\left(\phi \right)} \cos{\left(\theta \right)}\\- b \sqrt{1 - \frac{1}{R^{2}}} \sin{\left(\theta \right)}\end{matrix}\right]$




```python
v2
```




$\displaystyle \left[\begin{matrix}- a \sqrt{1 - \frac{1}{R^{2}}} \sin{\left(\phi \right)}\\a \sqrt{1 - \frac{1}{R^{2}}} \cos{\left(\phi \right)}\\0\end{matrix}\right]$




```python
import numpy as np
```


```python
def tangent_cone(a,b, P):
    """
    Возвращает три вектора `c,v1,v2`, определяющие эллипс, по точкам которого эллипсоид касается конус с вершиной `P`.
    
    Точки r_tan эллипса вычисляются так:
    ```
    c, v1, v2 = tangent_cone(a,b, P)
    t = np.linspace(0, 2*np.pi, 100)
    r_tan = c.reshape(3,1) + v1.reshape(3,1)*np.cos(t) + v2.reshape(3,1)*np.sin(t)
    ```
    В этом случае r_tan - матрица из 3 строк и 100 столбцов. Первая строка - координата `x` точек, 
    вторая строка - координата `y`, третья строка - координата `z`.
    """
    assert P.ndim == 1
    assert len(P) == 3
    
    x,y,z = P
    R = np.sqrt(x**2/a**2 + y**2/a**2 + z**2/b**2)
    cos_theta = (z/b)/R
    # theta меняется от 0 до pi, поэтому sin(theta) >= 0
    sin_theta = np.sqrt(1 - cos_theta**2)
    
    sin_phi = (y/a)/(sin_theta*R)
    cos_phi = (x/a)/(sin_theta*R)
    
    center = np.array([a*cos_phi*sin_theta/R, a*sin_phi*sin_theta/R, b*cos_theta/R])
    
    r = np.sqrt(1 - 1/R**2)
    v1 = np.array([a*r*cos_phi*cos_theta, a*r*sin_phi*cos_theta, -b*r*sin_theta])
    v2 = np.array([-a*r*sin_phi, a*r*cos_phi, 0])
    
    return center, v1, v2
```


```python
# Эллипсоид вытянут вдоль Z
a = 1
b = 4
```


```python
# Вершина конуса на бисектрисе первого квадранта
P = np.array([5,5,5])
```


```python
c, v1, v2 = tangent_cone(a,b, P)

```


```python
t = np.linspace(0, 2*np.pi, 100)
r_tan = c.reshape(3,1)+v1.reshape(3,1)*np.cos(t) + v2.reshape(3,1)*np.sin(t)
```


```python
el_X, el_Y = np.meshgrid(np.linspace(-a, a, 200), np.linspace(-a, a, 200))
```


```python
el_Z = b*np.sqrt(1 - el_X**2/a**2 - el_Y**2/a**2, )

```

    <ipython-input-38-80c6e10389cd>:1: RuntimeWarning: invalid value encountered in sqrt
      el_Z = b*np.sqrt(1 - el_X**2/a**2 - el_Y**2/a**2, )



```python
import matplotlib.pyplot as plt
```


```python
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(projection='3d')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.view_init(30, 120)

ax.plot_wireframe(el_X, el_Y, el_Z, cstride=10, rstride = 10, alpha=0.3)
ax.plot_wireframe(el_X, el_Y, -el_Z, cstride=10, rstride = 10, alpha=0.3)
ax.plot(r_tan[0], r_tan[1], r_tan[2], color="black")
ax.scatter(*P)

```


```python

```


```python

```

85b7053e4cfa71e61d959171beb3d6197026f6a7
ipynb
Jupyter Notebook
math/1293488/Untitled7.ipynb
pakuula/StackOverflow
4f3f37f98963980f592e2a1f5ee9165477166821
MIT
math/1293488/Untitled7.ipynb
pakuula/StackOverflow
4f3f37f98963980f592e2a1f5ee9165477166821
MIT
math/1293488/Untitled7.ipynb
pakuula/StackOverflow
4f3f37f98963980f592e2a1f5ee9165477166821
MIT
Qwen/Qwen-72B
1. YES
2. YES