# 由光源发起的遮盖计算

我们解决了照射计算的基本模型关系，并能够根据手电的位置指向，在地表求取光斑。但是，前文使用的是设置探针求取场强的点求取，对于绘制地表的等值线包络图、求取地表包线的具体解析情况，就不够用了。使用单点的方法计算量大，且步长不容易控制。

本文给出基于向量旋转与交汇的计算算法。

# 1. 向量参数方程与椭球交点

对一个高度为H的标准椭球，其方程为：

$$
\frac{x^2}{(a+H)^2}+\frac{y^2}{(a+H)^2}+\frac{z^2}{(b+H)^2}=1
$$

为了便于展开，设 $$$G=\frac{1}{(a+H)^2}$$$，$$$K=\frac{1}{(b+H)^2}$$$

假设ECEF一个空间向量为

$$
\vec T={A+Bt, C+Dt, E+Ft}
$$
大写的字母均为常数，t为参数方程的参数，则可以使用Octave的符号功能求取t与椭球的交点:

```octave
clc
clear
pkg load symbolic
G = sym('G');
K = sym('K');
A = sym('A');
B = sym('B');
C = sym('C');
D = sym('D');
E = sym('E');
F = sym('F');
t = sym('t');
T=solve(G*(A+B*t)^2 + G*(C+D*t)^2 + K*(E+F*t)^2==1,t);
latex(sym(T))
```
获得两个解为：

$$
\left[\begin{matrix}- \frac{A B G + C D G + E F K}{B^{2} G + D^{2} G + F^{2} K} - \frac{\sqrt{- A^{2} D^{2} G^{2} - A^{2} F^{2} G K + 2 A B C D G^{2} + 2 A B E F G K - B^{2} C^{2} G^{2} - B^{2} E^{2} G K + B^{2} G - C^{2} F^{2} G K + 2 C D E F G K - D^{2}
E^{2} G K + D^{2} G + F^{2} K}}{B^{2} G + D^{2} G + F^{2} K}\\- \frac{A B G + C D G + E F K}{B^{2} G + D^{2} G + F^{2} K} + \frac{\sqrt{- A^{2} D^{2} G^{2} - A^{2} F^{2} G K + 2 A B C D G^{2} + 2 A B E F G K - B^{2} C^{2} G^{2} - B^{2} E^{2} G K + B^{2} G- C^{2} F^{2} G K + 2 C D E F G K - D^{2} E^{2} G K + D^{2} G + F^{2} K}}{B^{2} G + D^{2} G + F^{2} K}\end{matrix}\right]
$$

# 2. 向量的绕轴旋转

向量的绕轴旋转在绘制地表区域时非常有用。考虑很多方向图都是对称的，或者想控制投影仪上的画笔在地表绘制线条，则可以控制一个 HPR向量在 $$$F(\alpha ,\beta)$$$ 底片上进行游动，其转换到ECEF下与椭球的交点就是投影的位置。



使用这样的方法，我们可以快速生成复杂的图形。
![2](../images/example03.png)

罗德里格旋转公式（Rodrigues' rotation formula）：

$$\vec v'=\vec v \cos(sita) + (1-\cos(sita))\cdot(\vec k\odot \vec v)\odot \vec k + \vec k \otimes \vec v \cdot sin(sita);$$
其中，v 是被旋转的向量， k 是转轴。

相关函数：

```cpp
/*!
 * \brief vec_rotation 罗德里格旋转公式（Rodrigues' rotation formula）
 * \param v		旋转前向量
 * \param k		旋转轴
 * \param sita  转角
 * \param res	旋转后向量
 * \param rad 弧度true/度 false
 */
void vec_rotation(const double v[3],const double k[3],const double dsita, double res[3],const bool rad = false);

/*!
 * \brief ellipsoid_intersections 求取过A的三维斜率为k的直线与椭球的交点。
 * \param A_ecef  A点的ECEF坐标
 * \param k  直线的单位矢量方向
 * \param H  椭球海拔
 * \param result  交点
 * \param pt  截距，即 A + kt 的 t
 * \return 交点个数。
 */
int ellipsoid_intersections(
	const double A_ecef[3],
	const double k[3],
	const double H,
	double result[/*2*/][3],
	double *pt/*2*/ = nullptr);
```