在上一篇中，我们解决了照射计算的基本模型关系，并能够根据手电的位置指向，在地表求取光斑。但是，前文使用的是设置探针求取场强的点求取，对于绘制地表的等值线包络图、求取地表包线的具体解析情况，就不够用了。使用单点的方法计算量大，且步长不容易控制。
本文给出基于向量旋转与交汇的计算算法。代码参考 https://gitcode.net/coloreaglestdio/geocalc
1. 向量参数方程与椭球交点
对一个高度为H的标准椭球，其方程为：
x 2 ( a + H ) 2 + y 2 ( a + H ) 2 + z 2 ( b + H ) 2 = 1 \frac{x^2}{(a+H)^2}+\frac{y^2}{(a+H)^2}+\frac{z^2}{(b+H)^2}=1 (a+H)2x2+(a+H)2y2+(b+H)2z2=1
为了便于展开，设 G = 1 ( a + H ) 2 G=\frac{1}{(a+H)^2} G=(a+H)21， K = 1 ( b + H ) 2 K=\frac{1}{(b+H)^2} K=(b+H)21
假设ECEF一个空间向量为
T ⃗ = A + B t , C + D t , E + F t \vec T={A+Bt, C+Dt, E+Ft} T=A+Bt,C+Dt,E+Ft
大写的字母均为常数，t为参数方程的参数，则可以使用Octave的符号功能求取t与椭球的交点:
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))
获得两个解为：
[ − A B G + C D G + E F K B 2 G + D 2 G + F 2 K − − 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 − A B G + C D G + E F K B 2 G + D 2 G + F 2 K + − 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 ] \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] [−B2G+D2G+F2KABG+CDG+EFK−B2G+D2G+F2K−A2D2G2−A2F2GK+2ABCDG2+2ABEFGK−B2C2G2−B2E2GK+B2G−C2F2GK+2CDEFGK−D2E2GK+D2G+F2K−B2G+D2G+F2KABG+CDG+EFK+B2G+D2G+F2K−A2D2G2−A2F2GK+2ABCDG2+2ABEFGK−B2C2G2−B2E2GK+B2G−C2F2GK+2CDEFGK−D2E2GK+D2G+F2K]
2. 向量的绕轴旋转
向量的绕轴旋转在绘制地表区域时非常有用。考虑很多方向图都是对称的，或者想控制投影仪上的画笔在地表绘制线条，则可以控制一个 HPR向量在 F ( α , β ) F(\alpha ,\beta) F(α,β) 底片上进行游动，其转换到ECEF下与椭球的交点就是投影的位置。
使用这样的方法，我们可以快速生成复杂的图形。
罗德里格旋转公式（Rodrigues’ rotation formula）：
v ⃗ ′ = v ⃗ cos ( s i t a ) + ( 1 − cos ( s i t a ) ) ⋅ ( k ⃗ ⊙ v ⃗ ) ⊙ k ⃗ + k ⃗ ⊗ v ⃗ ⋅ s i n ( s i t a ) ; \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′=vcos(sita)+(1−cos(sita))⋅(k⊙v)⊙k+k⊗v⋅sin(sita);
其中，v 是被旋转的向量， k 是转轴。
相关函数：
/*!* \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);
